Mon. Not. R. Astron. Soc. 000, 000-000 (1994) 



Printed 1 February 2008 



(MN plain TgX macros vl.6) 



On a matrix method for the study of small perturbations 
in galaxies 

Peter O. Vandervoort 

Department of Astronomy and Astrophysics, The University of Chicago, 5640 Ellis Avenue, Chicago, IL 60631, USA 

ABSTRACT 

A matrix method is formulated in a Lagrangian representation for the solution of the 
characteristic value problem governing modes of oscillation and instability in a colli- 
sionless stellar system. The underlying perturbation equations govern the Lagrangian 
displacement of a star in the six-dimensional phase space. This matrix method has a 
basis in a variational principle. The method is developed in detail for radial oscilla- 
tions of a spherical system. The basis vectors required for the representation of the 
Lagrangian displacement in this case are derived from solutions of the Lagrangian 
perturbation equations for radial perturbations in a homogeneous sphere. The basis 
vectors are made divergence free in the six-dimensional phase space in accordance 
with the requirement of Liouvillc's theorem that the flow of the system in the phase 
space must be incompressible. The basis vectors are made orthogonal with respect to 
a properly chosen set of adjoint vectors with the aid of a Gram-Schmidt procedure. 
Some basis vectors are null vectors in the sense that their inner products with their 
own adjoint vectors vanish. The characteristic frequencies of the lowest radial modes 
are calculated in several approximations for members of a family of spherical models 
which span a wide range of central concentrations. The present formulation of the 
matrix method can be generalized for nonradial modes in spherical systems and for 
modes in axisymmetric systems. 
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1 INTRODUCTION 

Considerations of stability and instability impose certain constraints on the equilibrium and evolution of a galaxy. Most 
galaxies observed in nature are presumed to be in stable states of equilibrium, or nearly so, inasmuch as they would not be 
expected to survive in unstable states long enough to be observed. If the processes of formation or the subsequent processes 
of evolution bring a galaxy into an unstable state, then the onset and development of that instability would act as dynamical 
mechanisms of evolution. Consequences of the instability could be manifest in a later equilibrium state of such a galaxy. 

Examples of instabilities in galaxies are well known. According to Toomre (1964), an axisymmetric disk would be locally 
unstable with respect to axisymmetric perturbations unless the dispersion of the peculiar velocities of the stars exceeded 
a certain critical value. The conjecture of Ostriker & Peebles (1973) that excessive rotation would make an axisymmetric 
galaxy unstable with respect to a deformation into a bar has become a central organizing principle in efforts to understand 
the dynamics of spiral galaxies. A third example, originally suggested by Antonov (1973) and demonstrated for the first time 
by Polyachenko & Shukhman (1981) and Polyachenko (1981), is the radial orbit instability in which the anisotropy of the 
velocity distribution makes a spherically symmetric stellar system unstable with respect to an axisymmetric deformation. 
Finally, certain firehose instabilities can deform the equatorial plane of an axisymmetric stellar system if the system is 
sufficiently flattened, and those instabilities might, according to Fridman & Polyacheko (1984; see Chapter X), explain the 
absence of elliptical galaxies flatter than type E7 on the Hubble sequence. 

It is beyond the scope of this introduction to present an extensive review of the history of investigations of such instabilities. 
An introduction to the subject is contained in the textbook of Binney & Tremaine (1987), and the subject is treated extensively 
in the monographs of Fridman & Polyachenko (1984) and Palmer (1994). Reviews by Merritt (1987, 1990) and Polyachenko 
(1987) are useful supplements to the monographs. Accounts of subsequent investigations (e.g., Merritt & Sellwood 1994; 
Sellwood & Merritt 1994; Robijn 1995; Sellwood & Valluri 1997) include references to more recent developments. 

For the work described in this paper, the most relevant part of the literature on small perturbations in galaxies concerns 
the matrix method formulated by Kalnajs (1977) for the solution of the characteristic value problem governing normal modes 
of oscillation and instability in a collisionless stellar system. Originally formulated for the study of instabilities in axisymmetric 
disks, that method has subsequently been used in several investigations of the radial orbit instability in spherical systems 
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(Polyachenko & Shukhman 1981; Fridman & Polyachenko 1984; Palmer & Papaloizou 1987; Weinberg 1989,1991a; Saha 1991, 
Bertin et al. 1994). More recently, the method has been extended to the investigation of small perturbations in axisymmetric 
systems (Robijn 1995). The matrix method of Kalnajs is based on an Eulerian representation of the perturbations. 

A Lagrangian representation of the perturbations provides an alternative framework for the investigation of the oscillations 
and the stability of stellar systems (Vandervoort 1983, 1989, 1991). This paper describes a matrix method for the study of 
normal modes in the Lagrangian representation. The elements of Lagrangian perturbation theory that provide the foundation 
for the method are reviewed in Section 2, and the method is formulated in Section 3. For the sake of definiteness and simplicity 
in this first presentation of the new matrix method, we concentrate on an application of the method to the study of radial 
modes of oscillation in a spherically symmetric galaxy. In Section 4 a procedure is described for the construction of basis 
vectors which can be used for the representation of such perturbations. The reduction and solution of the characteristic value 
problem in the matrix representation is described in Section 5, and Section 6 contains examples of solutions for modes which 
illustrate and test the method. A connection between the discrete and continuous spectra of modes in a stellar system is 
described in Section 7. Important technical details are described in a series of appendices. 

As this paper was being written, substantial progress was made in the generalization of the present formulation of the 
matrix method for a wider range of applications. The most important details of such generalizations have already been 
worked out for non-radial modes of oscillation in spherical systems and for general oscillations of axisymmetric systems. The 
construction of basis vectors for the representation of non-radial modes in spheres follows the procedure described in Section 4 
very closely. It consequently appears that much of the analysis described in this paper for radial oscillations in spheres can be 
readily extended to the study of non-radial oscillations in spheres. A similar construction of basis vectors for the representation 
of perturbations in axisymmetric systems can be achieved with the aid of a formulation based on an earlier investigation of 
the lowest modes of oscillation of a homogeneous spheroid (Vandervoort 1991). Finally, preliminary investigation suggests 
that this approach to the construction of basis vectors can be extended to nonrotating triaxial systems. The author plans 
next to apply these generalizations, in a continuation of the present work, to the radial orbit instability in spherical systems 
and to firehose instabilities in rotating axisymmetric systems. 



2 LAGRANGIAN PERTURBATION THEORY 

Small perturbations of a collisionless stellar system are described in the Lagrangian representation in terms of the Lagrangian 
displacement (Ax, Av) of a star (Vandervoort 1983). Thus, if the position and velocity of the star at time t would have been 
x and v, respectively, in the unperturbed system, then the position and velocity of the star at that time are, by definition, 
x + Ax and v + Av, respectively, in the perturbed system. 

For an infinitesimal perturbation, Aa; and Av are considered to be functions of x, v, and t, and the equations of motion 
for the perturbation are 

dAcc dAv . CA 1 ,-.s 

= At; and — — = Aa{Ax}, (1) 

where 

d _ 8 8 dV 8 

dt dt dx dx ' dv U 
is the total time derivative along the unperturbed trajectory of the star, and 

(3) 

n 

is the Lagrangian perturbation of the acceleration of the star. Here the distribution function fo(x,v) represents the unper- 
turbed density of stars in the six-dimensional phase space of a single star, and Vo(x) represents the gravitational potential in 
the unperturbed system. Also, G denotes the constant of gravitation, and m* denotes the mass of a single star. The second 
term on the right-hand side of equation (3) represents the Eulerian perturbation of the gravitational acceleration of a star, 
and the integration there extends over the region Q of the phase space that is accessible to stars in the unperturbed system. 

Equations (l)-(3) are a homogeneous system of linear integro-differential equations governing an infinitesimal perturba- 
tions of a stellar system. In Vandervoort (1989, hereafter LM), these equations have been reduced to a characteristic-value 
problem governing the normal modes of oscillation and instability of a system. The subsections that follow contain a review 
of the results derived in LM that are essential for the present formulation of a matrix method. 

2.1 The characteristic value problem 

For normal modes of oscillation having a time dependence of the form exp(— iu)t), where the characteristic frequency a; is a 
constant, we represent the Lagrangian displacement as a six-vector 



a r a \ ( a d\dV d d f Ax(x',v',t)f (x',v') , , 
Aa{Ax} = - (Ax — Gm,— - — dx dv 
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in the six-dimensional phase space, and we reduce equations (1) to the form 

lot] — Prj. (5) 
Here P is the matrix of operators 



-D i 
\iAa{} -D 



P 

where 

, , d dV 8 
D = i 



(6) 



/ d dV d\ 
V dx dx dv ) 

We identify r) as an element of a Hilbert space defined by letting the inner product of two six-vectors 
g (x,v)=( g f> V l) and Hx, V )=( h u ^ V \), (8) 



y 9 2 {x,v) J ' \h 2 (x,v) J ' 

say, where g x , g 2 , hi, and h 2 are three- vectors in the six-dimensional phase space, be 

(g,h} 2 =m* / g* (x, v) ■ h(x,v)f (x,v)dxdv, (9) 



where the asterisk, written as a superscript, signifies complex conjugation here and in what follows. In the notation adopted 
in LM, the subscript 2 appears on the left-hand side of equation (9) in order to distinguish the inner products of six-vectors 
from the inner products of three-vectors. The present definition of the inner product of two six-vectors differs slightly from 
the definition in LM; we now include a factor m, which was previously omitted. It is to be emphasized here that g* ■ h is the 
scalar product of two six-vector functions in the six-dimensional phase space. 

It is shown in the Appendix of LM that the Lagrangian displacement satisfies a divergence condition of the form 

V 6 ■ r) = A • Ax + A ■ Av = 0. (10) 
ox ov 

This condition implies, in accordance with Liouville's theorem, that the perturbed flow of the stars in the phase space is 
incompressible. 

2.2 The adjoint problem 

The operator P is not Hermitian in the Hilbert space of the Lagrangian displacements. Indeed, the operator adjoint to P is 
for it is shown in LM that, for any two six-vectors C,(x, v) and r](x, v), say, we have 

(C,Pr 1 ) 2 ^(P A C,v) 2 - (12) 

Inasmuch as the characteristic value problem represented by equation (5) is not self-adjoint, it must be solved in association 
with the solution of the adjoint problem 

* A nA A / i o\ 

u) r] — P r\ (13) 
governing vectors 



* A =W) (14) 

adjoint to the characteristic vectors rj. The notation on the right-hand side of equation (14) is intended to emphasize that we 
let the three-vectors Av A and Ax A have units of a velocity and a length, respectively, so that the unit of an inner product 
such as (C A 1*7)2' sa y> i s we ^ defined. 

As is explained in Section lid of LM, the association of the characteristic value problem with a distinct adjoint problem 
introduces some complication, in general, into the solution of the characteristic value problem. However, as is further explained 
in Section III of LM, the situation is greatly simplified in cases in which the structure of the unperturbed system is invariant 
with respect to the simultaneous application of time reversal and reflection through a plane. For the sake of definiteness we 
express that invariance by writing 

fo(xi,-x 2 ,x 3 , -vi,v 2 , -v 3 ) = fo(xi,x 2 ,x 3 ,v 1 ,v 2 ,v: i ), (15) 

where we are letting Xi and Vi (i = 1, 2, 3) denote the Cartesian components of x, and v, respectively, and we are letting the 
plane of reflection be the (xi, a;3)-plane. In this paper, we formulate and study the matrix method for systems having the 



Av' 
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symmetry described by equation (15). Examples of such systems include spherically symmetric systems, rotating axisymmetric 
systems in which the unperturbed distribution function fo{x, v) can be expressed in accordance with the theorem of Jeans as 
a function of the energy integral and the conserved component of the angular momentum of a star, and systems (including 
axisymmetric systems) with triplanar symmetry in which the gravitational potential Vo(x) is of the Stackel form and fo(x, v) 
is expressible in terms of three isolating integrals. 

The simplification of the characteristic value problem for such systems arises, because it is a consequence of equation(15) 
that a characteristic vector rf(u>) and its adjoint vector r] A (uj*) satisfy the relation 

V A (u*) = Ar,(u), (16) 



where A is the matrix 

C" 
C 



A 



(17) 



and the operator C transforms a three-vector function g(x,v,i) (say) in the phase space in the manner 

(gi(xi,X2,X 3 ,Vl,V2,V 3 ,t)\ / gi*(xi,-X2,x 3 ,-vi,V2,-v 3 ,-t) \ 

g2{xi,X2,x 3 ,vi,V2,v 3 ,t) -> Cg(x,v,t)= -g2*{xi,-x 2 ,x 3 ,-vi,V2,-v 3 ,-t) . (18) 
g 3 {xi,X2,X 3 ,Vl,V2,V 3 ,t) J \ g 3 *(xi,-X2,X 3 ,-Vl,V2,-V 3 ,-t) / 

The simultaneous application of time reversal, reflection through the (xi,x 3 )-p\axie, and complex conjugation transforms 
Ax and Av in the manner 

Ax(x, v, t) — > CAx(x, v, t) and Av(x, v, t) — ► — CAv(x, v, t), (19) 
respectively. The operator C commutes with the operators D and Aa{} in the sense that 

C{Dg) = D(Cg) and C(Aa{g}) = Aa{Cg}, (20) 
where g(x, v, t) is any three- vector function in the six-dimensional phase space. 

3 FORMULATION OF A RAYLEIGH-RITZ MATRIX METHOD 

We turn now to the matrix method for the solution of the characteristic value problem. The method is formulated here for 
normal modes whose frequencies form a discrete spectrum. 

3.1 A variational principle 

The basis for this matrix method is the following variational principle, which was proved in Section IIIc of LM. Let 



Aw J 

be an arbitrary six-vector in the six-dimensional phase space, and let the value of u> be determined by the relation 

(A V , (w - P) V ) 2 = 0- (22) 

Let 5ui be the change of u> that results when r] is subjected to an arbitrary variation Srj in equation (22). Then, under 
conditions in which (Ar),rj) 2 =fc 0, a necessary and sufficient condition that 5u> = is that ui and r) satisfy the equation 

lot) = P-q (23) 

and, accordingly, that r] is the characteristic vector belonging to the characteristic frequency lo. 

3.2 The matrix method 

Let the functions /3(x, v; k) (k = 1, 2, . . . .) be elements of a complete set of basis vectors in the Hilbert space of the Lagrangian 
displacements, and let these basis vectors satisfy the orthogonality conditions 

(Af3(j),(3(k))2 = if i^ k - (24) 

In view of equation (10), we require that the basis vectors also satisfy the divergence condition Ve • 0(k) = 0. We represent 
the Lagrangian displacement as a superposition of the basis vectors of the form 

»7 = ^a(fc)/3(fc), (25) 
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where the coefficients a(k) are constants. 

In order to apply the variational principle, we identify the right-hand side of equation (25) as a trial function for the 
Lagrangian displacement in which the coefficients a(k) are the variational parameters. Thus, we substitute from equation (25) 
for r\ in equation (22), and we make use of the definition of the inner product in equation (9), the definition of the operator 
A in equations (17) and (18), and the orthogonality relations in equation (24) in order to simplify the resulting equation. We 
next subject the resulting equation to mutually independent, arbitrary variations 5a(k) of the coefficients a(k). We require 
that the variation 8co of the frequency u> vanish in the variational equation that is obtained in this process. As a consequence 
of the independence and arbitrariness of the variations 5a(k), the variational equation can be satisfied, in general, only if the 
frequency u) and the coefficients a(k) satisfy the system of equations 



This is the matrix form of the characteristic value problem governing the normal modes of the system in the Lagrangian 
representation. Equations (26) are a homogeneous system of linear equations in the coefficients a(k). In general, a solution for 
those coefficients exists if and only if the determinant of the secular matrix of the system vanishes. That condition provides 
the characteristic equation for the determination of the frequencies ui of the modes. 

One might have thought that the matrix representation of the characteristic value problem could be derived directly from 
equation (5). However, it is not clear how, in such a direct approach, one would know that the orthogonality relations must be 
formulated in terms of inner products of the basis vectors and their adjoint vectors, that the matrix of the operator P is to be 
constructed as it is in terms of the basis vectors and their adjoint vectors, or that the matrix of the operator P is to be made 
symmetric as it is on the right-hand side of equation (26). On the other hand, these conditions are imposed automatically in 
the present derivation of the matrix equations from the variational principle. Thus the variational principle would appear to 
be an essential foundation for the matrix method. 

We have deliberately retained the inner products (A(3(k), f3(k)} 2 of the basis vectors and their adjoint vectors on the 
left-hand side of equation (26). Ordinarily, we might expect to normalize the basis vectors so that each of these inner products 
is equal to unity. However, we have found in practice, as is explained in Section 4 below, that the inner products of some 
basis vectors with their own adjoint vectors can vanish. In other words, (Af3(k), f3(j)} 2 is a diagonal matrix, in virtue of the 
orthogonality relations, with some vanishing elements on the principal diagonal. The non-vanishing elements on the principal 
diagonal may be normalized to unity in the usual way. 

Equations (26) are an infinite system of equations. In practice, we approximate the characteristic value problem by 
expressing the Lagrangian displacement as a superposition of a finite set of basis vectors /3(fc) on the right-hand side of 
equation (25) and truncating the system of equations (26) accordingly. In the following sections, we consider the solution of 
the characteristic value problem in such an approximation. 

4 A SET OF BASIS VECTORS FOR THE REPRESENTATION OF RADIAL OSCILLATIONS OF A 
SPHERICAL GALAXY 

In the remainder of this paper, we describe and illustrate the matrix method with the aid of an application to radial modes 
of oscillation of a spherical stellar system. Consistently with the theorem of Jeans, we let the distribution function of the 
unperturbed system be of the form fo(x, v) = fo(E,L 2 ), a function of the energy E and total angular momentum L of a 
star. Accordingly, the distribution function f (x,v) is an even function of the components of the velocity v. This section is 
devoted to a procedure for the construction of the basis vectors required for the representation of radial perturbations in such 
a sphere. 

The essential idea of the procedure is to construct a set of basis vectors which would suffice for the exact solution of the 
equations of motion governing the normal modes of oscillation of a particular model of a stellar system. The model chosen in 
the present case is a homogeneous sphere. Although the homogeneous sphere is an unrealistic model of a galaxy, the resulting 
basis vectors appear to give good results for the frequencies of radial modes even in quite centrally concentrated systems (see 
Section 6 below). 

For a homogeneous sphere of stars in which the unperturbed distribution function fo(x,v) — fo(E,L 2 ) is of a suitable 
form, the radial modes of oscillation have Eulerian perturbations of the gravitational potential which are polynomials in 
the square of the radial coordinate r — \x\ (see Chapter III, Section 4.1 of Fridman & Polyachenko 1984). The present 
construction of the basis vectors is accordingly based on solutions of the Lagrangian perturbation equations for the response 
of a homogeneous sphere to an imposed Eulerian perturbation of the gravitational potential which is a polynomial in r 2 . The 
solution for the Lagrangian displacement of the response is composed of six-vectors in the phase space which are polynomials 
in the Cartesian components of the position x and velocity v of a star. With a slight generalization for applications to systems 
with inhomogeneous density distributions, those 'six-vector polynomials' are almost the basis vectors that we are seeking. 
They are divergence free in the phase space, but they are not orthogonal in the sense of equation (24). The set of basis vectors 





(26) 
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finally adopted is constructed from the set of vector polynomials with the aid of a Gram-Schmidt procedure for satisfying 
the orthogonality conditions. Like the vector polynomials, the basis vectors are polynomials in the Cartesian components of 
x and v. 



4.1 Perturbations of a homogeneous sphere. 

For the purpose of investigating radial perturbations of a homogeneous sphere, it is convenient to reduce the Lagrangian 
perturbation equations as follows. In the interior of the sphere, the unperturbed gravitational potential is 

Vb(x) = f Gpor 2 = in V, (27) 
apart from an additive constant, where po is the density in the sphere, 

n 2 = ^-Gpo, (28) 

and r = \x\ is the radial coordinate of the point x in a system of spherical polar coordinates with the origin at the centre of 
the sphere. Stellar orbits in the unperturbed system are simple harmonic oscillations with frequency no. 

In the expression for Aa{Ax} in the second of equations (1) (see also eq. [3]), we substitute from equation (27) for Vo(x) 
and we identify 

Vl (x,t) = Gm*#- . / Ax(x',v' t)Mx',v') dx , dv , (2g) 



as the Eulerian perturbation of the gravitational potential (see eq.[14] in Vandervoort[1983]). We now eliminate Av between 
equations (1), and we bring the resulting equation of motion for Ax to the form 

In the integration of equation (30) in what follows, we shall make use of the requirement that x and v satisfy the unperturbed 
equations of motion 

dx i dv 2 , Q1 .. 

— = v and — = —no x. (31) 
at at 

For a radial oscillation of the system, the Eulerian perturbation of the gravitational potential is of the form Vi(x,t) — 
Vi(r) exp(— itdt), where u> is the frequency of the oscillation. When the radial function Vi(r) is imposed in the form of a 
polynomial in r 2 , the solution of equation (30) is a superposition of the solutions of a set of equations of the form 

d 2 \ 

— + no 2 J Ax = -xr 2k exp(-iwt) (k = 0, 1, 2 . . . .). (32) 

For the integration of equations (32), we introduce the variables 

z<j = x + aino~ v (a = +1, — 1). (33) 

After writing the right-hand side of equation (32) in terms of the new variables and expanding the resulting expression, we 
bring equation (32) to the form 

+ n„ 2 ) Ax = - (I) " £ £ £ 2- 1 f M f k ~ U \ z„M(z +1 , z_ i; fc - n - p, n, p) exp(-i^), (34) 

^ ' n=0 p=0 a \ J \ J 

where (™) denotes a binomial coefficient, the sum over a is over the values a = +1 and a = — 1, 

M(z+i,z-i;m,n,p) = (z+i ■ z +1 ) m (z +1 ■ z- 1 ) n (z- 1 ■ «-i) p , (35) 

and m(= k — n — p), n, and p are integers. Therefore, the solution of equation (32) can be written as a superposition of 
solutions of the equations 

d 2 2 \ 

^+n J Ax = -z a M(z + i,z-i;m,n,p)exp(-iu>t). (36) 
In virtue of equations (31), the variables z^ satisfy the equations of motion 
^- = -am z a (a = +1,-1). (37) 
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With the aid of equations (35) and (37), we readily verify that the particular solution of equation (36) for the response of the 
system to an imposed 'Eulerian perturbation of the acceleration' of the form — z a M(z+x, Z-i; m, n,p) exp(— iu)t) is 

z <7 M(z +1 ,z-x;m,n,p)exp(-iojt) 

Ax a (z + x,z-x,t;m,n,p) = — — — — ■ — r — ^ — . (38) 

[a; + (2m ~ 2p + a)n \ 2 - no 2 

The Lagrangian perturbation of the velocity is accordingly 

Av a (z+i, z-i,t; rn, n,p) = Ax„(z + x, Z-x, t; m,n,p) = -i[w + (2m - 2p + a)n ] Ax a (z +1 , Z-x, t; m, n,p). (39) 

The Lagrangian displacement described by equations (38) and (39) can now be represented as the six-vector 
' Ax a (z+i,z-i,t; m, n,p) 



(40) 



f Axa{z+i,z-i,t;m,n,p)\ 
V Av a (z+x, Z-x, t; m, n,p) J 

= 1 / z tr M(z+i,z-i;m,n,p) \ex (-iwt) 

[w + (2m — 2p + a)n ] 2 - n ' 2 V — i[w + (2m - 2p + o-)no]z a M(z+x, Z-x;m,n,p) ) 

For the purpose of writing the solutions of equations (32) , it is useful to rewrite equation (40) in the form 

CT/i+ 1 (« + i,z_i;m,n,p)exp(-io;t) a/^ 1 (z+i, z_i; m, n,p)exp(-iu;t) 

VcT{z+i,z-i,t;m, n,p) = — - — — — — — — — ■ 1 — ■ — — — —. - — ■ , (41J 

2no[u> + (2m — 2p + 2<j)no\ 2no\uo + (2m — 2p)no\ 

where the quantities 



(42) 



r, \ ( z a M(z+x,z-x;m,n,p) \ . . 

A* CT (*+i> 2 -i; m . n 'P) = • \ 1 (<t = +1,-1;t = +1,-1) 

yTcnnoZaM(z + i,z-i;m,n,p) J 

are 'six- vector monomials' in the three- vectors z a . Equation (41) represents the required solution of equation (36). Upon 
comparing equations (34) and (36), we see that the solution of equation (34) is 



{l\ k+1 S^S^S^nn ( k\ (k-n\ ( TCTnl(z +1 ,z- 1 ;k- rt - p, n,p) exp(-i^) \ 
V (z +1 ,Z- U t,k)--{-) ^^^2 III p II n()[u; + {2k _ 2n _ 4p + a + Ta)no] \, (43) 

n=0 p=0 a,r 



a superposition of the Lagrangian displacements represented by equation (41), where the sum over r is over the values r = +1 
and r = —1. In other words, equation (43) represents the Lagrangian displacement that satisfies equation (32) and describes 
the response of a homogeneous sphere to an 'Eulerian perturbation of the potential' equal to (2k + 2)~ 1 r 2fc+2 exp(— \ujt). 



4.2 Vector polynomials and their essential properties 

For k = 0, 1, . . ., we identify the vector polynomials to be used in the construction of the basis vectors by inspection of equation 
(43). For a given value of k, we observe that the possible values of the denominators on the right-hand side of equation (43) 
are lo+ (2k — 2n — 4p + o + To)no = lu + (2k + 2)no, cj + 2fcno, u> + (2k — 2)no, ■ ■ ., u> — (2fc + 2)no, apart from a common factor n . 
The number of distinct values of such denominators for a given value of k is equal to 2k + 3. Now we have obtained equation 
(43) as a solution of equation (30) (see also equation [32]) for an arbitrarily imposed Eulerian perturbation of the gravitational 
acceleration. In particular, the frequency to has an arbitrary value. On the other hand, the basis vectors and their ingredients 
should be independent of the value of u). Therefore, we suppress the factor exp(— iut) on the right-hand side of equation (43), 
and we identify the sum of the numerators there for each distinct value of the denominator ui + (2k — 2n — 4p + a + ra)no as 
a vector polynomial. There are thus 2k + 3 vector polynomials for a given power k of r 2 on the right-hand side of equation 
(32). This procedure determines the vector polynomials apart from arbitrary normalization factors. 
In the case that k — 0, for example, equation (43) can be written explicitly as 

77(2+1, t;0) 

1 / n+j(z + i,z- i; 0,0,0) nT_j(z +1 ,z_ i; 0,0,0) p±i(z + i,z-i;0,0,0) | fJ-lj (2+1, 0, 0, 0) \ ( 44 ) 
4no V w + 2no u> uj — 2no u> J ' 

where we have suppressed the factor exp(— iut) on the right-hand side. An inspection of the right-hand side of equation (44) 
as described in the preceding paragraph yields the three vector polynomials 

77(2+1, z-i; 0, 1) = z-i; 0,0,0), n(z+i, z_i; 0, 2) = fj+\{z+i, z-i; 0, 0, 0), and 

-1 -1 ( 45 ) 

77(2+1, 2_i ; 0, 3) = n_ 1 (z+i, z-x; 0, 0, 0) - n +1 (z+i, Z-x; 0, 0, 0). 

In the notation w(z+x,Z-x;k,h) for a vector polynomial, we adopt the convention that the integer k identifies the group of 
vector polynomials to which n(z+x, Z-x; k,h) belongs, i.e., the power k of r 2 on the right-hand side of equation (32), and 
the value of the integer h (h = 1, 2, . . . , 2k + 3) distinguishes 77(2+1, Z-x; k, h) from the other members of that group. In 
equations (A1)-(A3) in Appendix A, we list the expressions for the additional vector polynomials through the group k = 3. 
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We choose the normalization factors in the vector polynomials listed in equations (45) and (A1)-(A3) and we assign the order 
in which to list the vector polynomials in each group so that most of the vector polynomials are conveniently arranged in 
complex-conjugate pairs (see eq. [51] below). 

The vector polynomials listed in equations (45) and (A1)-(A3) would be ingredients for an exact construction of normal 
modes in the homogeneous spheres described by Fridman and Polyachenko (1984, see Chapter III, Section 4.1). In order to 
make use of the six- vector monomials /ij(m, n,p) and the six- vector polynomials Tr(k, h) in the representation of perturbations 
in inhomogeneous spheres, we must generalize the definition of no. A suitable generalization is provided by the equation 

no 2 = ~, (46) 
where 

W = ~J Po(r)r dV °^ dx and/ = J p (r)r 2 dx (47) 

denote the gravitational potential energy and the moment of inertia, respectively, of the unperturbed configuration, po(r) 
and Vo(r) are the density and gravitational potential, respectively, of the unperturbed configuration, and the integrations 
extend over the volume V of the unperturbed configuration. For homogeneous spheres, equations (28) and (46) are equivalent 
definitions of no. In general, however, the variables z a , the scalar monomials M(z+i, Z-i; m, n,p), the vector monomials 
/i£(z+i,z_i;m,n,p), and the vector polynomials ir(z+i,Z-i;k,h) (k = 0, . . . , 3) are defined by equations (33), (35), (42), 
and (45) and (A1)-(A3), respectively, with no defined by equation (46). 

Making use of equations (33), (35), and (42) in explicit calculations applied to the expressions for the vector polynomials 
listed in equations (45) and (A1)-(A3), we find that 

V 6 • n(k, h) = 0. (48) 

Thus, the construction of the basis vectors from the vector polynomials will ensure that the representation of the Lagrangian 
displacement in equation (25) will automatically satisfy equation (10). 

In the calculations that follow, we shall need the complex conjugates of the vector polynomials and the vectors adjoint 
to the vector polynomials. These are to be constructed in terms of the vector monomials, their complex conjugates and their 
adjoints. Referring to equations (33), (35), and (42), respectively, we find that 

z*=2- CT , M*(z + i,z-i\m,n,p) = M(z + i,z-r,p,n,m), and \ft T a {z+i,z-i;m,n,p)]* = ti T _ a (z+i,z-v,p,n,m), (49) 

where the asterisk, written as a superscript, labels the complex conjugate of a quantity. With the aid of equations (16)-(18), 
(33), and (49), we can also verify that C[z <r M(z+i, Z-i;m,n,p)] = z a M(z+i, Z-i; m, n,p) and, accordingly, that 

/0 C\ f z a M(z + i,z-i;m,n,p) \ / -Tain z a M(z + i, z-i; m, n,p) \ 
A„ a (z +1 ,z. 1 ;m,n,p)=^ c J y T<TinoZ<r M(z+i,z-i;m,n,p) ) = \ z a M(z +1 , z^m,n,p) | (50) 

The third of equations (49) implies that the vector polynomials listed in equations (45) and (A1)-(A3) satisfy the complex- 
conjugate relations 

n{k,h) = n*{k,h-l) {h = 2, 4, . . . , 2k + 2) and n(k, 2k + 3) = -w*{k, 2k + 3). (51) 

In other words, the group of vector polynomials for a given value of k contains k + 1 complex-conjugate pairs of members and 
one imaginary member. 

For the explicit construction of the basis vectors and for the associated numerical calculations, we relabel the vector 
polynomials in terms of a single integer. The prescription for doing so depends on the number of groups of vector polynomials 
to in used in the construction of the basis vectors. Thus, in the following subsection, we shall construct the basis vectors from 
the groups of vector polynomials n(k, h) for which k = 0, 1, . . . , fc max , say. 

We first relabel the vector polynomials which form complex-conjugate pairs by writing (see eqs. [51]) 

n(z + i, z-i; k, h) — 7r(z+i, Z-i; q), where q = k(k + 1) + h (k = 0, 1, . . . , Aw*; h = 1, 2, . . . , 2k + 2). (52) 

The value of the new integer q to be associated with a given pair of values (k, h) is determined here as follows. We recall that 
the group of vector polynomials 7r(/c', h'), for a given value of k' (k' < k), contains 2k' + 2 complex members. Summing that 
number from k' = to k' = k — 1, we count a total of k(k + 1) complex vector polynomials with < k' < k — 1. We then let 
q — k(k + 1) + h, inasmuch as the polynomial n(k,h) is number h in the listing of the polynomials belonging to the group 
k. The integer q takes the values q = 1, 2, . . . , Q in equations (52), where Q = (fc max + l)(fc m ax + 2) is the total number of 
complex vector polynomials in the groups such that < k < fc max . We next add the purely imaginary vector polynomials to 
the enumeration in equation (52) by writing 

n(z+i,z-i; k, 2k + 3) = n(z+i,Z-i; q), where q = Q + k + l (k = 0, 1, . . . , fc ma x). (53) 
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Here, the expression for g is determined by the circumstances that we are continuing the enumeration of vector polynomials 
begun in equation (52) and each group of vector polynomials contributes one purely imaginary member to the enumeration (see 
eq.[51]). The values of the integer q in equation (53) are q — Q + l,Q + 2, . . . ,N, where iV = Q + fc max + l = (fc ma x + l)(fcmax + 3) 
is the total number of vector polynomials, complex and purely imaginary, in the groups such that < k < fc max . 

Rewriting equations (51) in accordance with equations (52) and (53), we bring the complex-conjugate relations satisfied 
by the vector polynomials to the form 

7r(g) =7r*(g-l) (g = 2,4,...,Q) and n(q) = -n*(q) (q = Q + 1, Q + 2, . . . , TV). (54) 
According to equations (45) and (Al)-A3), the vector polynomials are superpositions of vector monomials of the form 



n(z+i,z-i;q) = ^(g|j)/4( 2; +i, 2; -i; m , n ,:P), 



(55) 



3 = 1 



where we label the vector polynomials 77(2+1, z-i; q) in terms of the single integer q in accordance with equations (52) and 
(53). On the right-hand sides of equations (55), we determine the constant coefficients (g|j), we determine the integers which 
distinguish the vector monomials n^(z+i, Z-i; m,n,p) as functions a — <?(q;j), t = r(q;j), m = m(q;j), n — n(q;j), and 
p — p(q;j), and we determine the number of terms J = J(q) in each sum by inspection of equations (45) and (Al)-A3). The 
coefficients (q\j) are seen there to be real. It follows that 



n*(z+i,z-i;q) = y^ j (q\j)[iJ. T a (z +1 , m, n,p)]* and An{z +1 , z_i; q) = ^^{q\j)Anl(z +1 , z_i; rn, n,p). 



(56) 



j=i 



4.3 Construction of the basis vectors 

We now write the basis vectors in the form 

"9-1 

0(1) = JV(1)tt(1) and (3{q) = N{q) 



^c(g,r)/3(r) +7r(g) 



(g = 2,3,...,iV), 



(57) 



where the constants N(q) and c(g, r) are to be determined with the aid of a Gram-Schmidt process. In that process, we require 
that the basis vectors be orthogonal in the sense of equation (24). Thus, for example, when we form the inner products of 
Af3(s) (s < q) with both sides of the second of equations (57) and apply the orthogonality conditions, we obtain 

N(q)[{Af3(s),f3( S )) 2 c(q,s) + {Af3(s),TT(q)) 2 ] =0 and (A/3(q), /3(g)> 2 = N{q){AP(q), n(q)} 2 

(g = 2,3,...,A> = l,2,...,g-l). 



(58) 



For the reduction of equations (58) to forms suitable for the systematic evaluation of N(q) and c(g, s), we need the 
equations 

-9-1 

J2c(l,r)Al3(r) + An(q) (q = 2, 3, . . . , N) (59) 



A/3(l) = N*(l)An(l) and Af3(q) = N*(q) 



for the vectors adjoint to the basis vectors /3(g). We obtain equations (59) by transforming equations (57) with the aid of 
equations (17) and (18). Equations (57) describe a representation of the basis vectors in which the constants N(q) must not 
vanish. Consequently the second of equations (59) allows us to reduce the first of equations (58) (except when q — 2 and 
s = 1) to 



(A0(s),l3(s)) 2 c(q,s) = -N( S ) 



c(«, r){Af3(r),n(q)) 2 + {An(s),n(q)), 



(g = 3,4,...,iV; S = 2,3,...,g-l). (60) 



The first of equations (58), again with the common factor N(q) suppressed, enables us to rewrite equation (60) as 



(Af3(s) J (3(s)) 2 c(q,s) = N(s) 



(A/3(r),l3(r)) 2 c(s, r)c(q, r) - (An(s), 7r(g)) s 



(g = 3,4,...,iV; S = 2,3,...,g-l). (61) 



Likewise, with the aid of the second of equations (59) and then the first of equations (58), we reduce the second of equations 
(58) to the form 

q-i 

(^7r(g),7r(g)) 2 -^(A/3(r),/3(r)) 2 c 2 (g,r) (q = 2, 3, . . . , N). (62) 



(Ap(q),P(q)) 2 = N 2 (q) 
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Finally, making use of the first of equations (59) in order to rewrite the first of equations (58) in the case that s = 1 and then 
the second of equations (58) in the case that q — 1, we obtain 

{Af3(l),f3(l)) 2 c(q,l) = -N(l){An(l),n(q)) 2 (q = 2, 3, . . . , TV) and (A/3(l), /3(1)> 2 = N 2 (1)(Att(1), tt(1)) 2 (63) 

as special cases of equations (61) and (62), respectively. The inner products (An(q'), 7r(g)) 2 in equations (61)-(63) are given 
quantities determined by the unperturbed structure of the system (see eq.[9]). Expressions required for the evaluation of 
(Air(q'),ir(q)) 2 are derived in Appendix B, and the symmetries and other important properties of (Air(q'), n(q)} 2 are described 
there. 

Equations (61)-(63) are the basic equations for the determination of the constants N(q) and c(g, s) in the representation 
of the basis vectors introduced in equations (57). Before we can solve those equations, we must normalize the basis vectors. 
In what follows, we shall find that the appropriate normalization is 

{A(3(q),(3(q)) 2 = 1 (g = l,2,...,Q) 

(64) 

= (q = Q + l,Q + 2,...,N). 

The unusual result that some of the inner products (Af3(q), /3(g)) 2 must vanish complicates the matrix method somewhat, 
but it seems not to create fundamental difficulties in the application of the matrix method. 

With the normalization of the basis vectors given, we can solve equations (61)-(63) for the quantities N(q) and c(g, s) in 
the following sequence. We first evaluate N(l) and c(g, 1) (g = 2, 3, . . . , N) in accordance with equations (63), and we next 
evaluate iV(2) in accordance with the first of equations (62). We subsequently solve equations (61) and (62) in the sequence 
g = 3, 4, . . . , N. For each value of g, in turn, we evaluate the constants c(g, s) in the sequence s = 2, 3, . . . , g — 1 with the aid 
of equations (61), and we finally evaluate N(q) in accordance with equation (62). At each step in this process, the values of 
c(s, r) and N(s) required for the evaluation of the right-hand side of equation (61) or the coefficient of N 2 (q) on the right-hand 
side of equation (62) are known as results of earlier steps in the process. 

In Sections 4.3.1 and 4.3.2 which follow, we present separate discussions of the solutions of equations (61)-(63) for the 
quantities N(q) and c(g, s) involved in the representations of the basis vectors satisfying the different normalization conditions 
in equations (64). The solutions for N(q) and c(g, s) satisfy certain identities which are derived in Appendix C. Some of the 
results described in Sections 4.3.1 and 4.3.2 and in Appendix C follow, because the solutions of equations (61)-(63) have the 
property that the quantities N 2 (q), c(g, s)/N(s), c(s,r)c(q,r), and c 2 (g, r) are all imaginary. For a proof of this result, we 
note that the inner products (Att(s), 7r(g)) 2 are imaginary quantities, as is shown in Appendix B, and we observe that we 
have normalized the basis vectors so that the inner products (Af3(s) , f3{s)) 2 are real quantities. Noting that equations (63) 
establish the results claimed here for N(l) and c(g, 1), we can prove the general result by induction with the aid of equations 
(61) and (62). 



4-3.1 Complex-conjugate pairs of basis vectors 



When < q < Q and < q < Q, the inner products {Ait(q'), 7r(g)) 2 do not vanish, in general. With the normalization given 
in equations (64), we find that we can satisfy equations (61)-(63) consistently when < q < Q. In particular, we can solve 
equations (63) in the manner 



iV 2 (l) = [<A7r(l),ir(l)) 2 ] 1 and c(g, 1) = -N(l)(An(l), n(q)} 2 (g = 2, 3, . . . , Q). 
We can also bring the subsets of equations (61) and (62) with q < Q to the forms 



c(g, s) =N(s) 
and 



J~] c(s, r)c(g, r) - (An{s),n(q))r 



(g = 3,4,...,Q;s = 2,3,. 



N\q) 



9-1 



(A7r(g),7r(g)} 2 -^c 2 (g,r) 



(g = 2,3,...,Q), 



(65) 



(66) 



(67) 



respectively. Equations (65)-(67) enable us to evaluate the constants c(g, s) and N(q), in the case that g < Q, in the sequence 
prescribed above in the paragraph following the normalization of the basis vectors in equations (64). 

An important consequence of equations (65)-(67) is that the basis vectors /3(g) (g < Q) form complex-conjugate pairs in 
the arrangement 



/3(g)=/3*(g-l) (g = 2, 4, . . . , Q). 



(68) 
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The expressions for the basis vectors given by equations (57) satisfy equations (68) because the vector polynomials Tt(q) satisfy 
equations (54) and the quantities determined by equations (65)-(67) satisfy the identities 



N(q) =N*(q-l) and c(q,q-l)=0 {q = 2, 4, . . . , Q) and 
c(q, s) = c*(q — 1, s — 1) and c(q, s — 1) = c*(q — 1, s) (q = 4, 6, 
These identities are verified in Appendix C. 



,Q;s = 2,4,. 



2). 



(69) 



4-3.2 Null basis vectors 

It remains to solve equations (61)-(63) in the case that Q < q < N. As an essential part of that solution, we must show that 
the resulting basis vectors are null vectors as we have claimed in equation (64). 

When s < Q < q < N, the normalization of the basis vectors (Af3(s), /3(s)) 2 = 1 (s = 1,2, ... ,Q) reduces the first of 
equations (63) and equations (61) to 



c{q,l) = -N(l)(An(l),n{q)) 2 
and 



(q = Q + l,Q + 2,...,N) 



c(q,s) = N(s) 



c(s,r)c(q,r) - (Air(s),ir(q))r 



(q = Q + l,Q + 2,...,N;s = 2,3,...,Q), 



(70) 



(71) 



respectively. Starting with the known values of c(q, r) and N(q), evaluated, as described Section 4.3.1, with the aid of equations 
(65)-(67), we can make use of equations (70) and (71) in order to evaluate the additional quantities c(q, s) in a sequence such 
that, at each step in the process, the values of quantities required for the evaluation of the right-hand side of equation (71) 
are known as results of preceding steps. 

We must finally consider the solution of equations (61) and (62) in the case that Q < s < q < N. In this case, the 
inner products (An(s), Tt(q)} 2 and {An(q),n(q)} 2 vanish according to equations (B7). The further reduction of equations (61) 
and (62) requires additional identities, especially equations (Cll), which are also derived in Appendix C. In virtue of the 
normalization (Af3(q), f3(q)) 2 = 1 (q = 1, 2, . . . , Q) and equations (Cll), the sums of terms from r = 1 to r = Q vanish on 
the right-hand sides of of equations (61) and (62). Thus, the subset of equations (61) and (62) which we must finally consider 
reduces to 

s-1 

(A(3(s), I3{s)) 2 c{q, s) = N(s) (Af3(r),f3(r)) 2 c(s, r)c{q, r) (q = Q + 2, Q + 3, . . . , N; s = Q + 1, Q + 2, . . . , q - 1), (72) 

r = Q + l 

{Af3(Q + l),f3(Q + l)) 2 =0, (73) 
and 

9-1 

{A(3(q),(3{q)) 2 = -N 2 (q) ^ (Af3(r), (3{r)) 2 c\q, r) (q = Q + 2, Q + 3, . . . , N). (74) 

r = Q + l 

According to equation (74), if (A/3(r), /3(r)) 2 = (r = Q + l, ...,q-l), then {A/3(q),/3(q)) 2 = 0. Equation (73) shows that 
the required condition is satisfied for q — Q + l. Therefore, by induction, we have (A(3(q), (3(q)} 2 = (q = Q + 1, Q + 2, , . . . , N) 
as is claimed in equation (64) . Consequently, equations (72) and (74) leave the values of the remaining quantities c(q, s) and 
N(q) indeterminate. Without loss of generality, we assign those values in the manner 



c(q, s) = {q = Q + 2, Q + 3, . . . , AT; s = Q + 1, Q + 2, . . . , q - 1) 
and 

N(q) = l + i (q = Q + l,Q + 2,...,N). 

The values of N(q) adopted in equations (76) preserve the property that the quantities N 2 (q) are imaginary. 



(75) 
(76) 



4.4 The completeness of the basis vectors 

The question arises as to whether or not the set of basis vectors constructed here is complete. We have not established that 
basis vectors constructed in this way form a complete set in general. We can only say that the adopted set of basis vectors 
is 'sufficiently complete' for the exact representation of radial modes in certain homogeneous spheres. Accordingly, one might 
hope that they would suffice for the representation of corresponding modes in inhomogeneous spheres. It appears in Section 
6 that they do indeed suffice for the representation of the lowest radial modes in rather centrally concentrated spheres. 
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The validity of the present formulation of the matrix method as a method for the approximate solution of the characteristic 
value problem governing modes does not require that the adopted basis vectors form a complete set. For the matrix method is 
a particular realization of a variational method in which the basis vectors are the ingredients of a trial function. The essential 
condition to be satisfied by the basis vectors is that they provide a sufficiently accurate representation of the Lagrangian 
displacements of the modes to be investigated. This is essentially the point of view adopted previously by Saha (1991). 

4.5 Construction of the matrix equations 

The construction of a suitable set of basis vectors reduces a given application of the matrix method to the evaluation of 
the matrix elements in a truncated set of equations (26) and then to the solution of the matrix equations. We describe the 
construction of the required matrix elements in Appendix D. 



5 SOLUTION OF THE CHARACTERISTIC VALUE PROBLEM IN THE MATRIX 
REPRESENTATION 

In the remainder of this paper, we shall restrict ourselves to the consideration of radial oscillations in systems in which the 
distribution function is a function of the energy alone and the velocity distribution is isotropic. The quadratures required for 
the evaluation of inner products and matrix elements are described in Appendix E. In what follows, we describe the reduction 
and solution of the matrix equations. 

5.1 Reduction of the matrix equations 

With the truncation of the set of basis vectors to a finite set containing JV members, and in virtue of the orthogonality and 
normalization of the basis vectors specified in equations (24) and (64), respectively, the matrix formulation of the characteristic 
value problem represented in equation (26) now reduces to 

N JV 

o;a(fc) = ^i? (Ar) (fc,j)a(j) (k = 1, 2, . . . , Q) and = ^ R {N) (k, j)a(j) (k = Q + 1, Q + 2, . . . , JV), (77) 

2 = 1 2 = 1 

where 

R (N) (k,j) = 1 -[(A(3{j),P(3{k))2 + (A0(k),PPU)U (78) 

We reduce the system of equations (77) to a standard form by making use of those members of the system in which k = 
Q + 1, Q + 2, . . . , TV in order to eliminate the amplitudes a(k) (k = Q + 1, Q + 2, . . . , JV) from the remaining equations in the 
system. We accomplish this with the aid of a sequence of transformations in each of which we eliminate one amplitude. After 
JV — K steps in that process, where K (Q + 1 < K < JV) is an integer, we will have reduced equations (77) to a system of the 
form 

K K 

ua(k) = ^R( K \k,j)a(j) (fc = l,2,...,Q) and = ^ R (K) (k, j)a(j) (k = Q + 1, Q + 2, . . . , K ), (79) 

2=1 2=1 

say, where the matrix is derived from the matrix R^ in a sequence of transformations of the form of equation (82) 

below. Equation (82) is the transformation with which we eliminate a(K). Thus, from the last of equations (79), we derive 

K-l 
2 = 1 

Upon substituting from this solution for a(K) in equations (79) (k = 1, 2, . . . , K — 1), we obtain the system 

K-l K-l 

o;a(fc) = ^i? (K - 1) (fc,j)a(j) (fc = 1, 2, . . . , Q) and = ^ R^^ (k, j)a(j) (k = Q + 1, Q + 2, . . . , K - 1), (81) 

2=1 2=1 

where 

**-»M-^M- ^<>^™ . (82) 

At the beginning of this sequence of transformations we have K = JV, and equations (79) coincide with equations (77). At 
the end of the sequence, we have K = Q + 1, and equations (81) reduce to 

Q 

o;a(fc) = ^i? w) (fc,j)a(j) (k = 1, 2, . . . , Q). (83) 

2 = 1 
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We are thus left with the characteristic value problem for the Q x Q matrix R^\k, j), where Q = (fc max + l)(fc max + 2). 

By definition (see eq. [78]) the matrix R^- N \k,j) has the symmetry R^- N \j, k) = R^ N \k,j). It follows that the matrices 
defined by equations (82) have the symmetry R {K \j,k) = R ( - K \k,j). In particular, R^ Q \j,k) = 7? (Q) (fc, j). Moreover, in 
the numerical investigation of the solutions of equation (83), we find that the matrix R^\k,j) has the following additional 
properties to the accuracy in which that matrix is evaluated. 

(i) For even values of k — j, the elements R^'(k,j) are real. 

(ii) For odd values of k — j, the elements R ( - Q \k,j) are imaginary. 

Thus the complex matrix R(Q\k,j) can be expressed in terms of a real matrix S^(k,j), say, in the manner 
R (Q \k,j)=i^S (Q \k,j). (84) 
If follows from this definition of S^\k,j) and the symmetry of 7?' Q '(fc, j) that 

S (Q) (k,j) = (-i) (k - j) R (Q) (k,j) and S- w) (j,fc) = (-l) (fe - 3) 5 (Q) (fc,j). (85) 
Letting 

a(k)=i h a(k) (k = 1,2,..., Q), (86) 
where the quantities a(k) are constants, we can now transform the system of equations (83) to the system 



va(k)=Y,S iQ) (k,j)a(j) (k = 1, 2, . . . , Q), (87) 



and we are left with a characteristic value problem for the real matrix S^\k,j). 

The numerical calculations required for the evaluation of the matrix S^^fc, j) and the solution of equations (87) have 
been performed with the aid of programs written in FORTRAN 77 and run in single-precision and complex arithmetic on 
Macintosh II and Macintosh Ilex computers. The program for the solution of equations (87) incorporates subroutines from 
the NAPAK library. Numerical quadratures required in the evaluation of inner products and matrix elements (see Appendix 
E) were performed with the aid of the fourth order Runge-Kutta integrator described by Press et al. (1986). 

These calculations have been performed in the four approximations defined by setting fc max equal to 0,1,2, and 3. The 
corresponding values of Q = (fc max + l)(fc max + 2) are 2, 6, 12, and 20, respectively. In other words, we have solved equations 
(87) for 2 by 2, 6 by 6, 12 by 12, and 20 by 20 matrices S^(k,j), respectively. For an assigned value of fc max , the number of 
basis vectors included in the representation of the Lagrangian displacement is TV = (fc max + 1 ) (fc max + 3) . Thus, we include 3,8, 
15, and 24 basis vectors in the four approximations considered here. The numerical calculations required for the evaluation of 

(k, j) and the solution of equations (87) in all four approximations require about four hours for one unperturbed model. 

5.2 Characteristic vectors 

Let (n = 1, 2, . . . , Q) denote the solutions of equations (87) for the characteristic frequencies for a given value of Q, and 

let a(j, n) (j — 1,2, ... ,Q;n = 1,2, ... ,Q) denote the corresponding solutions for the constants a(j). The solution for the 
Lagrangian displacement which belongs to the characteristic frequency Lo(n) is then 

JV 

ij(n)=5^a(j,n)/3CO (88) 
j=i 

(cf. eq. [25]). Here, 

a(j,n)=i j a(j,n) (j = 1,2,..., Q) (89) 

(cf. eq. [86]), and the remaining values of the constants a(j, n) (j = Q + 1, Q + 2, . . . , N) are determined by working backward 
(i.e., in the sequence K = Q + 1, Q + 2, . . . , N) through equations (80)-(82) with the aid of the values of the first Q constants 
a(j,n). 

Equation (88) represents a characteristic solution for the Lagrangian displacement in the approximation that the right- 
hand side of equation (25) is truncated to a superposition of TV basis vectors. These characteristic solutions satisfy orthogonality 
relations of the form 

(Ar)(m),r)(n)) 2 = 5 m „, (90) 

where <5 m „ is the Kronecker delta. In order to verify equations (90) , we first substitute from equations (88) for the Lagrangian 
displacements in the inner product on the left-hand side of equation (90) and bring the resulting expression to the form 

Q 

{Ar){m),n{n)) 2 = ^a{j,m)a{j,n) (91) 
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with the aid of equations (24) and (64). Now the requirement that the quantities ui(n) and a(j,n) (j = 1,2, ... ,Q;n — 
1,2, ... ,Q) be a characteristic solution of equations (87) implies that the quantities w(n) and a(j, n) (j — 1,2, ... ,Q;n — 
1,2, ... ,Q) must be a characteristic solution of equations (83). In the latter system of equations, the matrix R Q (k,j) is 
symmetric in the indices k and j. It therefore follows from equations (83) that the constants a(j, n) satisfy the conditions 

Q 

^ a(j, m)a(j, n) = 8 mn . (92) 

3 = 1 

Without loss of generality, we have imposed here the normalization condition that, when m = n, the sum on the left-hand 
side of equation (92) is equal to unity. Equations (91) and (92) combine to establish equations (90). These are precisely of 
the form of the orthogonality relations that must be satisfied by the exact solutions of the characteristic value problem (see 
eqs.[24] in LM). 

In virtue of equations (89), we can also combine equations (91) and (92) in the manner 

Q 

{Ar){m),ri{n)) 2 = ^(-l) J a(j, m)a(j, n) = S mn . (93) 

In other words, the constants a(j, n) are normalized in the manner 
Q 

^(-l)>(j,n)] 2 = l. (94) 

3 = 1 



5.3 Identification of modes in different approximations 

We shall distinguish solutions of equations (87) obtained in different approximations with the aid of subscripts. Thus, let 
ui a (n), a a (j,n), a a (j,n), and r) a {n) denote the quantities involved in the representation of the nth characteristic solution of 
equations (87) in the ath approximation, where a — fcmax + 1. Accordingly, we write equation (88) for the nth characteristic 
solution for the Lagrangian displacement in the ath approximation as 

N a 

V a (n)=Yla a (j,n)P a (j), (95) 

3 = 1 

where N a = (fc max + 1) (fcmax + 3) = a(a + 2). The basis vectors f3 a (j) are labelled here with a subscript a, because the sets of 
basis vectors constructed in Section 4.3 are different in the different approximations. The number of characteristic solutions 
of equation (87) in the ath approximation is equal to Q a = (fcmax + 1) (fcmax + 2) = a(a + 1) 

A useful comparison of the characteristic solutions r} a (n) and r) b {m) (a > b) is provided by the inner product 

N a N b 

Iab(n,m) = (Ar) a (n),ri b (m)} 2 = ^^a a {j,n)a b {k,m){Al3 a {j), (3 b (k)) 2 . (96) 

The values of \I ab (n, m)\ provide criteria for identifying characteristic solutions of equations (87) that represent the same mode 
in different approximations. Note that we take the modulus of the inner product here, because the phases of the solutions 
r} a (n) and r] b (rn) need not be the same. According to equation (24) in LM, the Lagrangian displacements of different modes are 
mutually orthogonal. If rj a (n) and r) b (m) were accurate representations of the Lagrangian displacement of the same mode, then 
we would have \I a b(n,m)\ = 1. Otherwise, if f] a (n) and ri b (m) were accurate representations of the Lagrangian displacements 
of different modes, then we would have \I ab (n,m) \ = 0. The characteristic solutions of equations (87) give only approximate 
representations of the Lagrangian displacements of modes, in general, so the values of \I ab (n,m)\ should lie between and 1 
in practice. 

Consider now the evaluation of the inner products (A/3 a (j), f3 b (k)) 2 on the right-hand side of equation (96). The pre- 
scription for the construction of basis vectors described in Section 4.3 (see particularly the ordering of the vector polynomials 
in equations [52] and [53] there) implies that 

f3 b (k) = f3 a (k) (fc = l,2,...,Q 6 ). (97) 
Therefore, we can reduce equation (96) to the form 

Qb N a N b 

I ab (n,m) = y^ a a (k,n)a b (k,m) + ^ ^ a a (j,n)a b (k,rn)(Af3 a (j), (3 b (k)} 2 (98) 

k=l j=Q b + lk=Q b + l 
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with the aid of the orthogonality relations governing the basis vectors (see eqs. [64]). Substituting from equations (57) for 
/3 b (k), and again making use of equations (64) and (97), we verify that 

{Af3 a (j),(3 b (k)) 2 = N b (k){Al3 a (j),n b (k)) 2 (Q b <j<N a ;Q b <k< N b ). (99) 
It also follows from equations (52) and (53) that 

TTb(k) = TT a (l) (k > Q b ;l = k + Q a - Q b >Q a ). (100) 
We also have 

N b (k) = N a {l) = l + i (k>Q b ;l = k + Q a -Q b >Q a ) (101) 

in virtue of equations (76). We can therefore reduce the right-hand side of equation (99) in the manner 

(Af3 a (j),(3 b (k)) 2 = N a (l){A/3 a (j),n a (l)) 2 

^ (102) 
= -N a (l) 22 c a (l,r){A/3 a (j),/3 a (r)) 2 (Q b < j < N a ; Q a < I = k + Q a ~ Q b < N b + Q a - Q b ), 

r=l 

where we have eliminated Tr a (l) with the aid of equations (57) and we have suppressed terms that vanish identically with 
the aid equations (24), (64), and (75). Upon substituting from equation (102) for (A(3 a (j), f3 b (k)) 2 on the right-hand side of 
equation (98) and making further use of equations (24) and (64) in the evaluation of the inner products (Af3 a (j), f3 a (r)} 2 , we 
finally obtain 

Qb Qa N b +Q a -Q b 

I ab (n,m) =*^2a a (k,n)a b (k,m) - ^ ^ a a (j,n)a b (l - Q a + Q b ,m)N a (l)c a (l, j). (103) 
fc=i i=0«+i 



6 NUMERICAL EXAMPLES 



In this section, we illustrate the matrix method with a calculation of modes of radial oscillation in a family of spherically 
symmetric models of galaxies. In each member of the family, the unperturbed distribution function is a function of the energy 
alone, and the unperturbed density is of the form 

3M r c 2 3M r c 2 



Po{r) = 



(r < R), 



(104) 



47r(r c 2 + r 2 ) 5 / 2 \_-k(t 2 + R 2 ) 5 / 2 

where R is the radius of the configuration, the constant r c is of the nature of a core radius, and the constant M is a certain 
characteristic mass. It follows from equation (104) that the characteristic mass is determined in terms of the central density 
po(0), the radius R, and the core radius r c in the manner 



Mo = ypo(0)# a :r c 3 



1 - 



(x 2 + 1)5/2 



where x c — r c /R. We can evaluate the mass of the system interior to a spherical surface of radius r as 

M r 3 M r 2 r 3 



M (r) = 4tt / p (r)r dr = 



( r . c 2 +r 2)3/2 ( rc 2 + #2)5/2 ' 



(105) 



(106) 



and we can accordingly evaluate the gravitational acceleration in the unperturbed system with the aid of Newton's theorem. 
In the limit R — » oo, the members of this family of models tend to Plummer's model. 

We make use of 7?, po(0), and [7rGpo(0)] _1/ ^ 2 as the units of length, density, and time, respectively. In this system of units, 
the models form a one-parameter sequence. The dimensionless core radius x c = r c /R is a convenient parameter with which 
to distinguish the members of the sequence. Alternatively, we can make use of the central concentration, defined here as the 
ratio C = po(0)/po of the central density to the mean density of a configuration, as the parameter which distinguishes models 
along the sequence. With the aid of equations (105) and (106), it is a straight-forward calculation to verify that 



C = 



Mo) 



YPo(0)i? 3 



4tt 



po(r)r dr 



1 



[(Xc 2 + if' ~ lei- 



(107) 



The solutions of equations (87) occur in pairs, the frequencies of each pair being equal in magnitude but opposite in 
algebraic signs, so it will suffice to present results of the present application of the matrix method only for modes with 
positive frequencies. Moreover, we consider only modes which we can identify in both the third and fourth approximations. 
As is explained at the end of Section 5, we make use of the values of the inner products (n,m) as criteria for the 
identification of such modes. For a given unperturbed system, we have 12 solutions of equations (87) in the third approximation 
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Table 1. Identification of modes of radial oscillation in a system in which C = 11.335 
(x c = 0.55). 



Mode 


n 


m 


I 


\hz{n,m)\ 


|/4 2 (n,0l 


|/ 3 2(m,0l 




0)3(771) 


w 2 (0 


Rl 


19 


9 


5 


0.9999 


0.9994 


0.9995 


0.6143 


0.6144 


0.6162 


R2 


20 


12 




0.9088 






0.8628 


1.1280 




R3 


18 


11 


6 


0.9243 


0.8448 


0.7429 


1.3000 


1.3510 


1.4180 


R4 


17 


10 




0.9938 






1.6360 


1.6460 




R6 


15 


8 


4 


0.8928 


0.5608 


0.6499 


2.4700 


2.6310 


2.3510 


R9 


4 


2 




0.8666 






4.1300 


4.2430 





(fcmai — 2,Q — 12), and 20 solutions in the fourth approximation (fe max = 3, Q = 20). Thus, we can expect to identify up to 
six modes with positive frequencies by comparing solutions of equations (87) in the third and fourth approximations. For this 
purpose we have evaluated \Inz(n, m)\ for all pairs of values (n, m) (n = 1,2, . . . , 20; m = 1, 2, . . . , 12), but we have arbitrarily 
selected for further investigation only those combinations in which 1/43(71, m)| > 0.4. Making similar use of the inner products 
\l42(n,m)\ and 1/32(71, m)\ as criteria, we can also identify solutions of equations (87) in the second approximation which 
represent up to three of the modes identified in the third and fourth approximations. 

Table 1 illustrates the identification of six modes in an unperturbed configuration in which x c = r c /R = 0.55 and 
C = po(0)/po = 11.335 (see eq. [107]). In the first column there, we label modes in terms of their characteristic frequencies in 
the fourth approximation. We arrange the positive solutions of equations (87) for the frequencies in the fourth approximation 
in a sequence in the order of increasing values, and we identify the nth frequency in that sequence and the mode to which it 
belongs with the label 'Rn'. The integers I, m, and n in the next three columns of Table 1 identify solutions of equations (87) 
in the second, third, and fourth approximations, respectively, with labels which were assigned by the LANPAK subroutines 
used in the solution of equations (87). The inner products 14,3(11, m), 142(11, 1), and 132(111,1) have been evaluated with the 
aid of equation (103). In the last three columns of the table, u>2(l), 0)3(771), and 0)4(71) denote solutions for the characteristic 
frequencies in the second, third, and fourth approximations, respectively. 

The solutions listed in Table 1 and their counterparts with negative frequencies are all of the cases in which |/4 3 (n, m)\ > 
0.4. The result that the listed values of I/43 (n, m) | (values which also apply to the corresponding cases with negative frequencies) 
exceed this limit by such a wide margin and approach unity so closely indicates that this comparison of the solutions of 
equations (87) in the third and fourth approximations is a reliable basis for the identification of modes. The trends in the 
values of \l32(m, l)\, \l42(n,l)\, and 1/43(71, m)\ for the modes which are also represented in the second approximation indicate 
that the sequence of approximations represented in Table 1 tends to converge. The values of the characteristic frequencies 
listed in Table 1 in the different approximations also indicate that the identification of modes in the different approximations 
is reliable and that the sequence of approximations tends to converge. 

For each of the modes represented in Table 1, we plot in Figures 1-6 the solutions in the second, third, and fourth 
approximations for the amplitudes a(k) and for the quantity 

' \x\ <r J \x\< 

which is the Eulerian perturbation of the mass interior to a spherical surface of radius r. In writing this expression for Mi(r, t), 
we have represented the Eulerian perturbation of the density pi(x, t) in terms of the hydrodynamical Lagrangian displacement 
€,(x,t) with the aid of the relations 

d f 

pi(x,t) = ■ (po£) and p (r)£(x, t) — m„ / Ax(x, v, t)f Q (x, v)dv (109) 



Mi(r,£) = Mi(r,0)exp(-iwt) = / pi(x,t)dx = - 

J\x\<r J\x\ 



Ax = -Atti- 2 p () (r)^ -£(x,t), (108) 



(Vandervoort 1983). In the second of equations (109), the integration extends over the region of the velocity space accessible 
to stars in the neighborhood of the point x. In order to evaluate £(x,t), we express the Lagrangian displacement rj as a 
superposition of the vector monomials (j,^(z + i, Z-\; m, n,p) with the aid of equations (55), (57), and (88), we extract from 
that superposition an expression for Ax(x,v,t), and we express the result of integrating Ax(x,v,t) over the velocity space 
as a superposition of the moments G(r, p) of the velocity distribution described in Appendix E. 

It happens that we can choose the phase of the solution for each mode so that the solutions for a(k) are real; we adopt 
that choice in all of the calculations described below. A comparison of the amplitudes a(k) in different approximations is a 
direct comparison of solutions of the characteristic value problem in the reduced matrix representation described by equations 
(87). Implicitly, this is a comparison of Lagrangian displacements in the different approximations. With the phase of the mode 
chosen so that the coefficients a(k) are real, it turns out that Mi(r, 0) is a complex quantity with equal real and imaginary 
parts. It will therefore suffice to plot the real part of M\(r, 0) in what follows. A comparison of the solutions for Re[Mi(r, 0)] 
in the different approximations is a comparison of the different approximations for the structures of perturbations in the 
configuration space. 

The mode Rl listed in Table 1 is the fundamental mode of radial oscillation. Evidently, the fundamental mode is repre- 
sented very accurately in these calculations. The plots of the solutions for a(k) and Re[Mi(r, 0)] are displayed in Figure 1. 
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Figure 1. Representations of the mode Rl in three approximations in a system in which C = 11.335 (x c = 0.55). This is the fundamental 
mode of radial oscillation, (a) Amplitudes a(k) in the representation of the perturbation, (b) Dependence of the Eulerian perturbation 
of the mass variable Rc[Mi(r, 0)] on the distance r from the centre of the unperturbed system. 



Figure 2. Representations of the mode R2 in two approximations in a system in which C = 11.335 (x c = 0.55). (a) Amplitudes a(k) in 
the representation of the perturbation, (b) Dependence of the Eulerian perturbation of the mass variable Rc[Mi(r, 0)] on the distance r 
from the centre of the unperturbed system. 

The mode is already well represented in the second approximation, and the additional amplitudes a(k) that appear in the 
third and fourth approximations are almost negligible. The solutions for Re[Mi(r, 0)] converge rapidly. 

The mode R2 is not represented very accurately in these calculations. According to Table 1, the discrepancy in the values 
of the frequency in the third and fourth approximations is substantial. In Figure 2a, the amplitudes a(k) (12 < k < 20) in 
the fourth approximation are not negligible, so the representations of the mode in the third and fourth approximations must 
be significantly different. The difference is manifest in the plots of Re[Mi(r, 0)] in Figure 2b. 

The convergence of the representations of mode R3 in the three approximations that is already evident in Table 1 is also 
evident in Figure 3. The amplitudes a(k) in the three approximations are in relatively good agreement, and, in the fourth 
approximation the amplitudes a(k) (12 < k < 20) are relatively small if not quite negligible. The plots of Re[Mi(r, 0)] show 
a satisfactory convergence. 

The representations of the mode R4 in Table 1 and Figure 4 appear to be quite accurate. The solutions for the amplitudes 
a(k) (1 < k < 12) in the third and fourth approximations agree very well, and the values of a(k) (12 < k < 20) in the fourth 
approximation are quite small. The runs of Re[Mi(r, 0)] in the two approximations are remarkably similar. 

The values of \l32(m, l)\ and |/42 (n, l)\ listed in Table 1 for the mode R6 indicated that that mode is not accurately 
represented in the second approximation. In Figure 5(a), the amplitudes a(k) in the second approximation do not agree very 
well with their counterparts in the third and fourth approximations, and the amplitudes a(k) (k > 6) are substantial in the 
third and fourth approximations. Likewise, in the runs of Re[Mi(r, 0)] in Figure 5(b), the second approximation does not 
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Figure 3. Representations of the mode R3 in three approximations in a system in which C = 11.335 (x c = 0.55). (a) Amplitudes a(k) 
in the representation of the perturbation, (b) Dependence of the Eulerian perturbation of the mass variable Re[Mi (r, 0)] on the distance 
r from the centre of the unperturbed system. 



Figure 4. Representations of the mode R4 in two approximations in a system in which C = 11.335 (x c = 0.55). (a) Amplitudes a(k) in 
the representation of the perturbation, (b) Dependence of the Eulerian perturbation of the mass variable Re[Mi(r, 0)] on the distance r 
from the centre of the unperturbed system. 



provide an accurate representation of the mode. Nevertheless, the agreement of the representations of the mode R6 in the 
third and fourth approximations is quite good. 

In the plots of a(k) and Re[Mi(r, 0)] in Figure 6, it appears that the identification of the mode R9 in Table 1 is reliable 
and that the representation of the mode is relatively accurate. 

The preceding discussion illustrates the matrix method and identifies six of the lowest modes of radial oscillation in 
a system with a relatively modest central concentration (C = po(0)/po = 11.335; x c = r c /R — 0.55). Figure 7 shows the 
dependence of the characteristic periods P — 2n/u> on the central concentration C for four modes which can be identified in 
both the third and fourth approximations and traced over a range of substantially greater central concentrations. Values of the 
periods derived in the second, third, and fourth approximations are plotted as open triangles, circles, and squares, respectively. 
We have chosen to plot the periods instead of the frequencies in Figure 7 in order to facilitate comparisons of the different 
approximations in systems with higher central concentrations. The criterion adopted here for the identification of a mode in 
the third and fourth approximations is the condition \Iiz{n,m)\ > 0.75. When this criterion is satisfied, points representing 
the periods in the third and fourth approximations are plotted as large circles and large squares respectively. The smaller 
circles and squares represent extrapolations or interpolations to values of the central concentrations at which the values of 
I/43 (n,m) | are less than 0.75. Two of the modes presented in Figure 7 are also represented in the second approximation, and 
those results are plotted as large triangles or small triangles according as |7 42 (n, Z)| > 0.75 or |J 42 (n,i)| < 0.75. 

In the present calculations, it has not been possible to obtain accurate solutions of equations (87) in the fourth ap- 
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Figure 5. Representations of the mode R6 in three approximations in a system in which C = 11.335 (x c = 0.55). (a) Amplitudes a(k) 
in the representation of the perturbation, (b) Dependence of the Eulerian perturbation of the mass variable Re[Mi (r, 0)] on the distance 
r from the centre of the unperturbed system. 



Figure 6. Representations of the mode R9 in two approximations in a system in which C = 11.335 (x c = 0.55). (a) Amplitudes a(k) in 
the representation of the perturbation, (b) Dependence of the Eulerian perturbation of the mass variable Re[Mi(r, 0)] on the distance r 
from the centre of the unperturbed system. 

proximation for unperturbed systems in which C > 1025 (x c < 0.1). The failure of these calculations for the most centrally 
concentrated systems considered was traced to the accumulation of round-off errors in the numerical solution of equations 
(61)-(63) for the quantities N(q) and c(q,r). 

The results displayed in Figure 7(a) indicate that the present representations of the fundamental mode, the mode Rl, 
remain quite accurate in systems with central concentrations up to values at least of the order of 10,000. The behavior of the 
characteristic frequency in the fourth approximation is described in the most centrally concentrated systems by the formula 

uj « 3.5(7rGp ) 1/2 (C* >> 100), (110) 

where pa is the mean density of the system. For the fundamental mode in systems in which C = 11.34 (x c = 0.55) and in 
which C — 1025 (x c = 0.1) plots of a(k) and Re[Mi(r, 0)] are shown in Figures 1 and 8, respectively. Whereas the second, 
third, and fourth approximations are substantially the same in the less centrally concentrated system represented in Figure 1, 
the three approximations are rather different in the more centrally concentrated system represented in Figure 8. Nevertheless, 
the sequence of approximations does appear to be converging to an accurate representation of the mode in the latter case. The 
node near the origin in the run of Re[Mi(r, 0)] in Figure 8(b) is unexpected. That feature need not be real; it could be a result 
of the accumulation of numerical errors mentioned above or it could be an artifact of the present method of approximation. 
That this feature is absent in the second approximation, present in the third approximation, and almost absent in the fourth 
approximation suggest that it might be an artifact. The feature does not appear at lower central concentrations. 
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Figure 7. Dependence of the characteristic period P = 2n/u> on the central concentration C for four modes of radial oscillation. Open 
triangles, circles, and squares represent the second, third, and fourth approximations, respectively. The different sizes of the symbols are 
explained in the text. Filled circles represent the results of N-body solutions of the Lagrangian perturbation equations. 

Figure 7(a) includes a comparison of the present results with values of the period of the fundamental mode calculated for 
the same unperturbed systems with the aid of an N-body method for the solution of the Lagrangian perturbation equations 
(Vandervoort 1999). The N-body results, which are represented by filled circles in Figure 7(a), were obtained before the start of 
the present work on the matrix method. A mode is excited in such N-body experiments by imposing suitable initial conditions 
on the Lagrangian displacements (Ax, Av) of the particles. In the experiments represented in Figure 7(a), suitable initial 
conditions were found by a process of trial and error in which corrections were applied to virial estimates of the Lagrangian 
displacements (Vandervoort 1983,1999). Those efforts to excite the fundamental mode without exciting other perturbations 
significantly were successful only for systems in which the central concentrations are sufficiently low (C < 137.8; x c > 0.2). 
In future work, it will be of interest to make use of Lagrangian displacements derived with the aid of the matrix method as 
initial conditions with which to excite the fundamental mode and other modes in N-body experiments on systems covering a 
wider range of central concentrations. 

In Figure 7(b), it appears that the identification of the mode R2 is quite firm over the full range of central concentrations 
spanned by the plotted points. Even the smaller points plotted at C = 17.20, 22.49, and 31.01 represent cases in which 
\U.i(n,l)\ > 0.7. Figure 9 shows plots of a(k) and Re[Mi(r,0)] for the mode R2 in the case that C = 1025 (x c = 0.1). The 
plots of a(k) in the third and fourth approximations appear to be quite different; nevertheless the runs of Re[Mi(r, 0)] in 
the two approximations agree rather well. The node that appears near the origin in the run of Re[Mi(r, 0)] in the fourth 
approximation seems to be an artifact of the approximation like the similar feature mentioned above in the case of the 
fundamental mode. The fourth approximation to the mode R2 does not have such a feature when C = 137.8 (x c = 0.2). A 
comparison of Figure 9(b) with Figure 2(b) indicates that the mode R2 has significantly different structures in systems of 
high and low central concentrations. The run of Re[Mi(r, 0)] has two interior nodes at C = 11.34 (x c = 0.55) but only one 
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Figure 8. Representations of the mode Rl in three approximations in a system in which C = 1025 (x c = 0.1). This is the fundamental 
mode of radial oscillation, (a) Amplitudes a(k) in the representation of the perturbation, (b) Dependence of the Eulerian perturbation 
of the mass variable Rc[Mi(r, 0)] on the distance r from the centre of the unperturbed system. 



Figure 9. Representations of the mode R2 in two approximations in a system in which C = 1025 (x c = 0.1). (a) Amplitudes a(k) in 
the representation of the perturbation, (b) Dependence of the Eulerian perturbation of the mass variable Re[A4"i(r, 0)] on the distance r 
from the centre of the unperturbed system. 

interior node at C = 1025 (x c = 0.1). The transition from two nodes to one node occurs at C « 13.59 (x c ~ 0.5) in the third 
approximation and at C « 17.20 (x c « 0.45) in the fourth approximation. 

Although it was shown above that, for the system in which C = 11.34 (x c — 0.55), the identification of the mode R4 is 
reliable and the representation of the mode is accurate, the identification and representation of the mode becomes problematic 
at higher central concentrations. In the representation of the mode in Figure 7(c), the runs of two of the solutions for the 
characteristic frequency in the third approximation are shown. The application of the criterion 1/43(71, m)| > 0.75 implies 
that the lower of the two curves represents the mode R4 in the cases that C < 22.49 (x c > 0.4) whereas the higher curve 
represents the mode in the cases that 31.01 < C < 137.8 (0.35 > x c > 0.2). At central concentrations exceeding 137.8, we find 
no solution of equation (87) in the third approximation that can be identified convincingly as a representation of the mode 
R4. (Of course the nomenclature adopted in order to label modes implies that the mode R4 is correctly represented in the 
fourth approximation by definition.) Plots of a(k) and Re[Mi(r, 0)] are shown in Figure 10 for the mode R4 in the case that 
C = 137.8 (x c = 0.2). These results in the third and fourth approximations agree well enough to support the identification of 
the mode in both approximations but not well enough to show that we have an accurate solution for the mode in this case. 
And we have no indication of the accuracy with which the mode is represented in the fourth approximation at greater central 
concentrations. 

As is shown in Figure 7(d), the mode R8 can be identified in accordance with the criterion |/43(n, m) > 0.75 only in 
centrally concentrated systems. The smaller circles and squares represent points in the third and fourth approximations, 
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Figure 10. Representations of the mode R2 in two approximations in a system in which C = 137.8 (x c = 0.2). (a) Amplitudes a(k) in 
the representation of the perturbation, (b) Dependence of the Eulcrian perturbation of the mass variable Rc[Mi(r, 0)] on the distance r 
from the centre of the unperturbed system. 



Figure 11. Representations of the mode R8 in two approximations in a system in which C = 1025 (x c = 0.1). (a) Amplitudes a(k) in 
the representation of the perturbation, (b) Dependence of the Eulcrian perturbation of the mass variable Re[Mi(r, 0)] on the distance r 
from the centre of the unperturbed system. 

respectively, at which 0.55 < \Ii3{n, m)\ < 0.75. Likewise, the small triangles represent points on the run of the period in 
the second approximation at which 0.55 < I/42 (w, Z) < 0.75. Plots of a(k) and Re[Mi(r, 0)] are shown in Figure 11 for the 
mode R8 in the case that C = 1025 (x c — 0.1). The solutions for the amplitudes a(k) in the three approximations agree 
quite well in Figure 11(a), but the amplitudes that are added at higher values of k in the third and fourth approximations are 
significant. The three approximations for the run of Re[Mi (r, 0)] agree quite well in Figure 11(b). It should be noted, however, 
that the three approximations do not agree on the number of nodes in the outer regions of the system. The hydrodynamical 
Lagrangian displacement is badly determined in the outer regions of a centrally concentrated system because there is very 
little mass there. In particular, the (two or) three approximations often disagree on the number of interior nodes in the run 
of the hydrodynamical Lagrangian displacement near the outer boundary of the unperturbed system. 

7 THE CONTINUOUS SPECTRUM OF NORMAL MODES 

In all of the cases described in Section 6, the modes that have been identified have frequencies which lie within the continuous 
spectrum of the frequencies of the radial components of the stellar orbits. Specifically, the frequencies of the modes considered 
here lie within the interval 



GMp(R) 
R 3 



, 2 - 
< uj < 



16-ivGpo 
3 



(111) 
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The lower limit in this inequality is the Keplerian frequency of an orbit at the outer edge of the system, whereas the upper 
limit is the frequency of the radial component of a (small) simple-harmonic oscillation of a star about the centre of the 
system. Inequality (111) implies that the modes identified here lie within the continuous spectrum of radial modes that is to 
be expected in the systems considered. 

In general, the normal modes of oscillation of a galaxy include a set of modes with a discrete spectrum of real and 
complex frequencies and a set of modes with a continuous spectrum of real frequencies. Modes of instability would belong to 
the discrete spectrum. The frequencies belonging to the continuous spectrum satisfy resonance conditions with the frequencies 
of the unperturbed stellar orbits in the system. In particular, frequencies satisfying inequality (111) would satisfy appropriate 
resonance conditions for a continuous spectrum of radial modes in the systems considered in the present work. Examples of 
systems in which the normal modes include both discrete and continuous spectra of frequencies are described by Fridman & 
Polyachenko (1984; see particularly Chapter II, Section 7 and Chapter III, Sections 3.4.1 & 6.1.3). Mathur (1990; see also 
Weinberg 1991b) has investigated conditions under which the normal modes of a system include a discrete spectrum of real 
frequencies which lie outside the continuous spectrum. 

The question arises in the present work as to how discrete modes can appear within the continuous spectrum of modes. 
The probable explanation is that such discrete modes are isolated members of the continuous spectrum with frequencies which 
'accidently' satisfy the characteristic equation. 

This interpretation is supported by an investigation of the solutions of the Lagrangian perturbation equations for the 
complete spectra of normal modes in stellar systems (Vandervoort 2001). That investigation is based on a matrix method for 
the study of perturbations which is rather different from the method described in this paper but similar to the matrix method 
of Kalnajs (1977). For our present purposes, the principal results are the following. For the discrete spectrum of modes in 
a given system, there is a characteristic equation for the determination of the frequencies. In general, however, there is no 
characteristic equation governing frequencies in the continuous spectrum. On the other hand, there are exceptional frequencies 
in the continuous spectrum which do satisfy the characteristic equation for discrete modes. Although the details of that work 
must be left to a later publication, the circumstances in which isolated frequencies in the continuous spectrum can satisfy the 
characteristic equation can be illustrated with the aid of well known results in the theory of plasma oscillations. 

In the second matrix formulation of Vandervoort (2001) , the modes of oscillation of a galaxy are recognized as counterparts 
of the van Kampen modes of oscillation in a homogeneous plasma (van Kampen 1955, 1957; Case 1959). For the sake 
of simplicity, we consider a model of a plasma which consists of an electron gas of uniform density in the presence of a 
fixed, uniform background of positive charge. A van Kampen mode is an electrostatic plane wave of wave number k and 
frequency to, say. The simplest example is obtained for a one-dimensional representation in which the electrostatic field is 
parallel to the direction of the wave, and functions of the velocity of an electron are integrated over the components of the 
velocity perpendicular to the direction of the wave. The Eulerian perturbation of the distribution function of the electrons 
and the electrostatic field are determined by solving Vlasov's equation (i.e., the linearized, collisionless Boltzmann equation) 
and Poisson's equation simultaneously. The Eulerian perturbation of the distribution function is of the form fi(x,v,t) — 
g(v)exp[i(kx — cot)] where the a;— direction is the direction of the wave and the amplitude g(v) is a function of the a;— component 
v of the velocity. The solution for g(v) is of the form 

g(v) = up 2 P^^-— — + XS(u-v), (112) 
dv u — v 

where fo{v) is the distribution function of the unperturbed electron gas, u = ui/k is the phase velocity of the wave, up = u}p/k 
is a characteristic velocity derived from the plasma frequency lu p , the symbol P denotes the Cauchy principal value, 6(u — v) is 
Dirac's delta function, and A is a constant to be determined. Applying the normalization condition J g{v)dv = 1 to equation 
(112), we obtain 



a=i-«p 2 p r *m 1 d „. aw) 

/ dv u — v 

J — c- 



For an arbitrary choice of the value of u = ui/k, in general, equation (113) determines the value of A. Thus, in general, there 
is no dispersion relation (characteristic equation) for the determination of the frequencies to of the van Kampen modes in the 
continuous spectrum. For isolated frequencies, however, it can happen that A = 0. Such frequencies can be determined by 
setting the right-hand side of equation (113) equal to zero. The condition that the right-hand side of equation (113) vanish 
is the dispersion relation originally derived by Vlasov (1945) for plasma oscillations. It has been remarked (Bohm & Gross 
1949, Bernstein, Greene & Kruskal 1957, Jackson 1960) that Vlasov's dispersion relation picks out the frequencies of waves 
in which there are no electrons trapped in the electrostatic field. 

In precisely the same way, isolated members of the continuous spectrum of modes in a galaxy are found to satisfy the 
characteristic equation for discrete modes. It is commonly expected that such perturbations would suffer Landau damping. 
However, Landau damping should not occur in individual modes in a spherical system in which the unperturbed distribution 
function depends only on the energy of a star. According to a theorem of Antonov (1961; Fridman & Polyachenko 1984; 
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Binney & Tremaine 1987; Palmer 1994), such a system will be stable if the distribution function is a monotonically decreasing 
function of the energy, as it would be in the systems considered in Section 6. Moreover, a lemma proved in LM shows that 
the frequencies of normal modes in a stellar system occur, in general, in complex-conjugate pairs. Thus, the frequencies of the 
modes considered in this paper must be real. In general, Landau damping does not occur in modes of the Van Kampen type. 
On the other hand, perturbations which are superpositions of modes belonging to the continuous spectrum would be expected 
to suffer Landau damping (Van Kampen 1955, 1957; Case 1959; see also Appendix 5. A. in Binney & Tremaine 1987). Our 
concentration on a normal mode analysis puts the study of Landau damping beyond the scope of the present investigation. 
The present calculations give no evidence for the existence in the systems studied of a discrete spectrum of radial modes with 
frequencies lying outside the continuous spectrum. 

In the N-body experiments on the fundamental mode of radial oscillation that are represented in Figure 7(a), the orbits 
that are populated with stars have frequencies of the radial motions which do not span the entire range described by inequality 
(111). These experiments were performed on systems of 1000 bodies. The density of stars is so low in the regions of the phase 
space accessible to orbits with the lowest frequencies that these N-body realizations of the systems studied do not include such 
stars. In other words, the N-body realizations of these systems do not include stars that would resonate with the fundamental 
mode. Thus the fundamental mode appears to be a discrete mode in the N-body experiments, and that perturbation should 
not suffer Landau damping. For substantially larger (and more realistic) numbers of bodies, N-body realizations would include 
stars that resonate with the frequency of the fundamental mode, and perturbations that include the fundamental mode could 
suffer Landau damping. The damping would be expected to be slow, because the density of resonant stars in the phase space 
would be very low. 

8 CONCLUDING REMARKS 

For the study of small perturbations in stellar systems, the matrix method described in this paper has two particularly 
attractive features. The first is that the characteristic value problem for modes in the matrix representation is linear in the 
frequencies of the modes. Consequently, the solution of the characteristic equation for the frequencies is readily obtained with 
the aid of standard matrix methods. In contrast, the matrix method of Kalnajs involves a characteristic equation which is 
nonlinear in the frequencies, and it is a more challenging task to obtain the solutions for the frequencies. The second feature 
that makes the present matrix method particularly attractive is that the underlying Lagrangian representation deals directly 
with the perturbations of the trajectories of the particles. The comparison of analytic investigations of small perturbations with 
N-body experiments is straight forward in this case. Analytic solutions for the Lagrangian displacements of the perturbations 
can be used as initial conditions for N-body experiments which provide direct tests of analytic results. The experimental 
results included in Figure 7(a) represent such a test. 

The results obtained in Section 6 consist of relatively accurate representations of six modes at central concentrations 
of the order of 10, quite an accurate representation of the fundamental mode of radial oscillation at central concentrations 
ranging up to the order of 10,000, and an approximate delineation of three other modes at central concentrations ranging 
up to the order of 100-1,000. These results were obtained with the aid of a very modest computational effort. Therefore, 
considerable scope remains for more diverse investigations of modes of oscillation and instability in stellar systems with the 
aid of more extensive numerical calculations. 

For these reasons, the planned extension of the present investigation to nonradial modes in spheres and to modes in 
axisymmetric systems, as described at the end of Section 1, would appear to hold considerable promise. 
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APPENDIX A: LIST OF VECTOR POLYNOMIALS 

In this appendix, we list additional groups of vector polynomials identified with the aid of the procedure described at the 
beginning of Section 4.2. 

The vector polynomials in the case that k = 1 are 

tt(1, 1) = 0, 0), tt(1, 2) = M ±i(0, 0, 1), 

tt(1, 3) = /4i(0, 1, 0) + |/*Zi(l, 0, 0) - 0, 0), tt(1, 4) = M ti(0, 1, 0) + ±n+{{0, 0, 1) - \liZ\(0, 0, 1), and (Al) 

tt(1, 5) = /4i(o, 0, 1) - 0, 0) + 2 M I 1 (0, 1, 0) - 2 M ;J (0, 1, 0), 

where we suppress the arguments z + i and z_i here and in what follows. 
Likewise, the vector polynomials in the case that k = 2 are 

tt(2, 1) = /n+i(2, 0, 0), tt(2, 2) = M li(0, 0, 2), 

tt(2, 3) = 1, 0) + \nZ\(2, 0, 0) - \n+\(2, 0, 0), tt(2, 4) = »±\(0, 1, 1) + ±it+{(0, 0, 2) - ^(O, 0, 2), 

tt(2, 5) = 0, 1) + 2 M +i(0, 2, 0) + 2nZ\{l, 1, 0) - \»Z\{2, 0, 0) - 2/^(1, 1, 0), (A2) 

tt(2, 6) = nZ\{\, 0, 1) + 2/*li(0, 2, 0) + 2 M ; 1 (0, 1, 1) - \n%\{Q, 0, 2) - 2/^(0, 1, 1), and 

7r(2, 7) = n+\(Q, 1, 1) + \ttZ\(l, 0, 1) + /xli (0, 2, 0) - n±\(l, 1,0)- \n+\^, 0, 1) - H+{(0, 2, 0). 
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Finally, the vector polynomials in the case that k = 3 are 
tt(3, 1) = M ^(3, 0, 0), tt(3, 2) = M ti(0, 0, 3), 

tt(3, 3) = n+\{2, 1, 0) + |/*Zi(3, 0, 0) - |M^i(3, 0, 0), tt(3, 4) = M ±i(0, 1, 2) + |^+^(0, 0, 3) - gfili (°, 0, 3), 
tt(3, 5) = /*+J(2, 0, 1) + 4/*+i(l, 2, 0) + 2 M l{(2, 1, 0) - |#*±i(3, 0, 0) - 2/^(2, 1, 0), 
tt(3, 6) = 0, 2) + (0, 2, 1) + 2 M ;J (0, 1, 2) - ifi+J(0, 0, 3) - 2/^(0, 1, 2), 

tt(3, 7) = 1, 1) + |/i+i(0, 3, 0) + |/*=i(2, 0, 1) + pC\{l, 2, 0) - ^11(2, 1, 0) - \/t+l{2, 0, 1) - ^(1, 2, 0), (A3) 

7r(3, 8) = M ^(l, 1, 1) + f Mti(0, 3, 0) + 0, 2) + ^(0, 2, 1) - ±/i+i(0, 1, 2) - 0, 2) - M :J(0, 2, 1), 

and 

tt(3, 9) = 0, 2) + 4 M +i(0, 2, 1) + 4^(1, 1, 1) + |m!i(0, 3, 0) 

- jitl (2, 0, 1) - 4 M li(l, 2, 0) - 4 M ; 1 1 (1, 1, 1) - |m;!(0, 3, 0). 

APPENDIX B: INNER PRODUCTS OF VECTOR POLYNOMIALS AND THEIR ADJOINTS 

The construction of the basis vectors described in Section 4.3 requires the evaluation of the inner products {An(q'),n(q)} 2 . 
These can be expressed in terms of the inner products (Afi^/(m', n',p'), V> T a {m, n,p))2- 

Making use of equations (42) and (50) for the vector monomials and their adjoint vectors, respectively, constructing the 
inner products of those vectors in accordance with equation (9), and taking the complex conjugates of quantities in accordance 
with equations (49), we obtain 

{Afi^,(m ,n ,p), fil(m,n,p))2 = ino(rV + To)m* I z_ iy , ■ z a M(m + p , n + n ,p + m')f dxdv. (Bl) 

Jn 

By inspection of equation (Bl), we see that 

(Afj^,{m ,n ,p'),fj^(m,n,p)}2 = unless to 1 = to. (B2) 

We also note that the quantity z_ CT / • z a M(m + p' , n + n ,p + m') can be expressed as a polynomial in the components of 
x and iv in virtue of equations (33) and (35). For systems of the kind considered in Section 4, in which fo(x,v) is an even 
function of v, it follows that the imaginary part of the integrand on the right-hand side of equation (Bl) is an odd function of 
the components of v. Thus the integral on the right-hand side of equation (Bl) is a real quantity, and we have the result that 
n',p'), Hv(m, n,p))2 is an imaginary quantity for systems in which the distribution function fo is an even function 
of the components of the velocity. Finally, we verify that 

{Avl(m,n,p),tf r ,(m',ri,p'))2 = (A^,(m',n',p'),^(m,n,p)> 2 . (B3) 

For, as we have just explained, the integral on the right-hand side of equation (Bl) is real. Therefore, we may replace the 
integrand there with its complex conjugate, rewrite the result with the aid of the first two of equations (49) and, in virtue of 
equation (Bl), transform an expression for the right-hand side of equation (B3) into an expression for the left-hand side of 
equation (B3). 

It follows from equations (55) and (56) that 

j' J 

(An(q'),n(q)) 2 = ]T ]>>'| j') (q\j){A^J, (m', n',p'), ^(m, n,p)) 2 , (B4) 

j' = l j = l 

where the values J = J(q), o = o(q;j), t = r(q;j), m = m(q-j) n = n(q;j), p = p(q;j), J' = J(q'), o' = o(q';f), 
t' — T(q';j'), m' — m(q';j') n' = n(q';j'), and p' = p(q';j') are determined by inspection of equations (45) and (A1)-(A3), 
as we have explained following equation (55). 

The inner products {Air(q'),n(q)) 2 have the following properties. 

(i) Inasmuch as the coefficients (q'\j') and (q\j) are real constants in equation (B4) and the inner products 
(Ap,^,(m' ,n' ,p'), nl(m,n,p)}2 are imaginary quantities, the inner products (Aw(q'),n(q)) 2 are imaginary quantities. 

(ii) Rewriting equation (B4) with the aid of equation (B3), we readily verify that {An(q'),n(q)) 2 has the symmetry 

(An(q'),n(q)) 2 = (An(q),n(q')) 2 (q' = 1, 2, . . . , N; q = 1, 2, . . . , N). (B5) 
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(iii) A consequence of the first of equations (54) is that 

(An(q'),n(q)) 2 = (An(q' - l),n(q - 1))* 2 and (An(q'),n(q - 1)> 2 = (An(q' - l),n(q))* 2 

(B6) 

(g' = 2,4,...,Q;g = 2,4,...,Q), 

inasmuch as complex conjugation commutes with the operator A. 

(iv) Likewise, it follows from the second of equations (54) that 

(An(q'),n(q)) 2 = (q' = Q + 1, Q + 2, . . . , N; q = Q + 1, Q + 2, . . . , N). (B7) 

More specifically, the second of equations (54) implies that the inner product on the left-hand side of equation (B7) is equal 
to its complex conjugate and is therefore real. However, the general result described in item (i) above implies that this inner 
product is imaginary. Therefore, the inner product described in equation (B7) must vanish. Equation (B7) is essentially a set 
of orthogonality relations for the purely imaginary vector polynomials and their adjoint vectors. In the case q' = q, we see 
that the inner product of a purely imaginary vector polynomial and its own adjoint vector vanishes. In other words, the purely 
imaginary vector polynomials are null vectors. It is for this reason that some of the basis vectors constructed in Section 4.3 
turn out to be null vectors as well. 

(v) A further consequence of equations (54) is that 

(An(q'), iz{q)) 2 = (An(q'), 7r(g - 1)> 2 {q = Q + 1,Q + 2, ... ,N;q = 2,4, ... ,Q). (B8) 

For, when we substitute from equations (54) for the vector polynomials on the left-hand side of equation (B8) and recall that 
complex conjugation and A are commuting operators and that these inner products are imaginary quantities, we obtain the 
inner product on the right-hand side of equation (B8). 

(vi) Finally, we observe that 

{An{q),n(q-l)) 2 = {An{q-\),n(q)) 2 = (q = 2, 4, . . . , Q). (B9) 

The first of the equalities in equation (B9) is a special case of equation (B5). When we combine that symmetry relation with 
the result of putting q — q in the second of equations (B6) , we find that the inner products in equation (B9) must be real. On 
the other hand, such inner products have been shown to be imaginary, in general, so they must vanish as claimed in equation 
(B9). Equation (B9) is essentially a set of orthogonality relations for the complex-conjugate pairs of vector polynomials and 
their adjoint vectors. 

Our procedure for the evaluation of the integral over the phase space on the right-hand side of equation (Bl) is described 
in Appendix E below. The numerical values of the inner products (ATr(q'), 7r(g)) 2 derived with the aid of equation (B4) have 
been checked with the aid of equations (B5)-(B9). 



APPENDIX C: CERTAIN IDENTITIES 

The arrangement of the basis vectors /3(g) (0 < q < Q) in complex-conjugate pairs as described in Section 4.3.1 (see eq. [68]) 
is a consequence of the identities presented in equation (69) . In what follows, we shall prove equations (69) by induction. 
We begin the proof by reducing the second of equations (65) and equation (67) in the case that q — 2 to 

c(2, 1) = -N(l){An(l), tt(2)> 2 = and iV 2 (2) = [(An(2), n(2)) 2 - C 2 (2, 1)] - 1 = [(An(l), ^(l)};]- 1 = [iV 2 (l)]*, (CI) 

respectively, with the aid of equations (54) and (B9) and the first of equations (65). It follows from the first of equations (CI) 
that we can reduce equation (66) in the case that s = 2 to 

c(q, 2) = N(2) [c(2, l)c(q, 1) - (An(2), n(q)) 2 ] = -N(2){An(2), 7r(g)} 2 (q = 3, 4, . . . , Q). (C2) 
Equations (B6) and the second of equations (CI) (see also eqs. [65]) enable us to derive the relations 
c( 9 ,2) = -JV*(l)(A7r(l),7r( 9 -l)> 2 *=c*(g-l,l) and c(q, 1) = -N* (2){An(2), n(q - 1)) 2 * = c (q - 1, 2) 

(g = 4,6,...,Q) 

from equation (C2). Equations (CI) establish the first two of equations (69) in the case that g = 2, and equations (C3) 
establish the third and fourth of equations (69) in the case that q = 4, 6, . . . , Q and s = 2. 

Continuing the proof of equations (69), we now consider equation (66) for given, even values of q and s, we recall that 
the vector polynomials satisfy equations (B6), and we assume, in the sense of a proof by induction, that N(s), c(s,r), and 
c(g, r) (r < s < g) satisfy equations (69). Under such conditions, we can rewrite equation (66) in the manner 



c(q,s) = N*(s-l) 



s-2 



^c*(s - l,r)c*(g- l,r) - (Att(s - l),7r(g - 1))* 



(C4) 



= c*(g - 1, s - 1) (g = 6, 8, . . . , Q; s = 4, 6, . . . , q - 2). 
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In this reduction, we break the sum on the right-hand side of equation (66) into sums over even and odd values of r, rewrite 
the separate sums with the aid of the third and fourth of equations (69), suppress c(q,q — 1) in accordance with the second 
of equations (69), and recombine the results in the single sum shown in equation (C4). In the verification of this derivation, 
it is helpful to note that the application of the third and fourth of equations (69) turns sums over even and odd values of r 
into sums over odd and even values of r — 1, respectively. Under the same conditions and with similar reductions, we can also 
rewrite equation (66) in the manner 



c(q,s-l) = N*(s) 



y^c*(s,r)c*(q - l,r) - (An(s),n(q - 1))* 



(C5) 



= c\q - 1, a) {q = 6, 8, . . . , Q; s = 4, 6, . . . , q - 2) 
and, in the case that s = q — 1, in the manner 

9-2 

c(q, q-l)=N{q-l)J2 <1 ~ 1. r)c(q, r ) 

7-1 

= N(q - 1) (q, r)c(q - 1, r) = (q = 4, 6, . . . , Q). 

r=l 

The quantity c(q, q — 1) must vanish here, because, the equality of the two sums in equation (C6) implies, on the one hand 
that each sum is a real quantity, whereas our discussion just prior to Section 4.3.1 shows, on the other hand, that the terms 
c(q — l,r)c(q,r) and c*(q,r)c*(q — 1, r) in those sums must be imaginary quantities. Finally, we can transform the right-hand 
side of equation (67) with the aid of equations (B6) and (69) and derive 

N 2 (q)= ({An(q-l),n( q -l)): 2 - q ^r[c 2 ( q -l,r)A = [N 2 (q - 1)}* (g = 4, 6, . . . , Q). (C7) 

The completion of the proof of equations (69) now consists of starting with the results in equations (CI) and (C3), 
considering equations (C4)-(C7) in the sequence q = 3, 4, . . . , Q, and, for each value of q, verifying equations (C4) and (C5) 
in the sequence s = 2, 3, . . . , q — 2 and then verifying equations (C6) and (C7). This confirms the arrangement of the basis 
vectors /3(g) (0 < q < Q) in complex-conjugate pairs as described in equation (68). 

The construction of the null basis vectors in Section 4.3.2 requires the reduction of a subset of equations (61) and (62) 
to equations (72)- (74). We conclude this appendix with the derivation of identities that are required for that reduction. We 
must first prove that the quantities c(g, s) which satisfy equations (70) and (71) also satisfy the identities 

c(q,s) = -c*(q,s-l) (q = Q + 1, Q + 2, . . . , N; s = 2, 4, . . . , Q). (C8) 

For a proof of equation (C8) by induction, we first consider equation (71) in the case that s = 2, and we make use of 
equations (54), (CI), and (70) in order to write 

c(q, 2) = -N(2)(An(2), n(q)) 2 = N*(l)(An(l), ir(q)) 2 =-c*(q,l) (q = Q + 1, Q + 2, . . . , AT), (C9) 

which proves the identity in the case that s = 2. In order to complete the proof, we must show that the condition c(q, r) — 
— c*(g, r — 1) holds for r — s, where s(> 2) is an even integer, if that condition holds for even values of r which are less than 
s. For this purpose, we make use of equations (54) and (69) in order to reduce equation (71) to 



c(q,s) = -N*(s-l) 



s-2 

^c*(s - l,r)c*(q,r) - {Air(s - l),7r(g)>* 

r=l 



(CIO) 



= -c*{q,s-l) (q = Q + 1, Q + 2, . . . , N; s = 4, 6, . . . , Q). 

In this reduction, we break the sum on the right-hand side of equation (71) into sums over even and odd values of r, rewrite 
the separate sums with the aid of the third and fourth of equations (69) and the condition c(q,r) = —c*(q,r — 1) (r = 
2, 4, . . . , s — 2), make use of the second of equations (69) in order to suppress c(s, s — 1), and recombine the results in the 
single sum shown in equation (CIO). That the transformed expression for c(q, s) in equation (CIO) is equal to — c*(q, s — 1) is 
also a consequence of equation (71). 

The identity that is required for the derivation of equations (72)-(74) is 

Q Q 

J2<s,r)c(q,r) = £V (a, r)c*(g, r) =0 (q = Q + 2, Q + 3, . . . , N; a = Q + 1, Q + 2, . . . , q). (Cll) 

r— 1 r— 1 

In this reduction, we break the first sum in equation (Cll) into sums over even and odd values of r, rewrite the separate 
sums with the aid of equations (C8), and recombine the results in the second sum shown in equation (Cll). The two sums 
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in equation (Cll) must vanish, because, their equality implies, on the one hand, that each sum is a real quantity, whereas 
our earlier discussion of equations (61)-(63) shows, on the other hand, that the terms c(s,r)c(q,r) and c* (s,r)c* (q,r) in those 
sums must be imaginary quantities. 



APPENDIX D: THE MATRIX OF THE OPERATOR P 



In order to evaluate the matrix (Af3(q'), Pf3(q)) 2 , we apply the operator P to the right-hand sides of equations (57), and we 
form the inner product of Af3(q') with the resulting expressions. We obtain 



(A/3(q'),P(3(l)) 2 = N(l)(A/3(q'),Pn(l)) 2 
and 



(A(3(q'),P(3(q)) 2 = N(q) 



9-1 



c(«, r){Af3(q'),Pf3(r)) 2 + {Af3(q'), Pn(q)) 2 



(<?>!)• 



(Dl) 



(D2) 



In order to evaluate the right-hand sides of equations (Dl) and (D2), we form the inner products of the adjoint vectors Af3(q'), 
expressed as in equations (59), with the quantities P-n(q). We obtain 



(49(1), Pn(q)} 2 = N(l){An(l), Pn(q)} 2 
and 

rv-i 

{A(3(q'),P7v(q)) 2 = N(q') 



c(q',r')(Af3(r'),Pn(q)) 2 + (An(q% Pn(q)) 2 



(q' > 1). 



Finally, we construct the matrix 

j' J 

(An(q'),Pn(q)) 2 = ^^^\j'){ q \j){A^l' l (m\n\p'),P^(m,n,p)) !l 



(D3) 



(D4) 



(D5) 



(see eqs. [55] and [B4]) in order to evaluate the right-hand sides of equations (D3) and (D4). 

Evidently, the calculations involved in the evaluation of (Aj3{q), P/3(l)} 2 begin with the evaluation of the matrix 
{A{i T a i(yri ,n' ,p'), Py, T a (m,n,p)) 2 and continue with the evaluation of (Aiv(q'), PTr(q)) 2 . We next consider equations (D3) 
and (D4) in the sequence q = 1, 2, ... , and, for each value of q in turn, evaluate the quantities {A/3(q'), Prt(q)} 2 in the se- 
quence q' = 1, 2, . . .. At each step, the right-hand side of equation (D3) or (D4) can be evaluated in terms of quantities whose 
values are known results of preceding steps in the procedure. Likewise, we finally consider equations (Dl) and (D2) in the 
sequence q — 1,2,..., and, for each value of q in turn, evaluate the quantities {Af3(q'), Pf3(q)) 2 in the sequence q' — 1,2, . . .. 
Similarly, the right-hand side of equation (Dl) or (D2) can be evaluated at each step in terms of quantities whose values are 
known results of preceding steps. 

We have now reduced numerical work required for the construction of the matrix (Af3{q ), Pf3(q)} 2 to the evaluation of 
the matrix 



{A(il,(m' ,n' ,p'),Pnl(m,n,p)) 2 — — ir' '& 'rcjn 2 m t I z_ a , ■ z a M(m + p ' ,n + n ' ,p + m')fodxdv 

Jn 

— i(r' g + rcr)nom,* / M(p , n , m')z_ cr r ■ D[zcrM(m, n,p)]fodxdv 
Jn 

+ im* / M(p\n\m')z_ a i ■ Aa{z a M(rn,n,p)} fodxdv. 
Jn 



(D6) 



In the derivation of equation (D6), we have constructed {Afi^t^m' ,n' ,p'), Pfi^(m,n,p)) 2 in accordance with the definition of 
an inner product given in equation (9), we have substituted from equations (50), (6), and (42) for A/x^, (jn! , n',p'), P, and 
li T a (m,n,p), respectively, in the integral to be evaluated, and we have reduced equation (D6) to the form given here with the 
aid of equations (35) and (49). 

The evaluation of the third integral on the right-hand side of equation (D6) requires the construction of the quantity 
Aa{zcrM(m,n,p)}. Let po(r), Mo(r), and Vo(r) denote the density, the mass interior to a spherical surface of radius r, and 
the gravitational potential, respectively in the unperturbed system. By definition and by Newton's theorem, we have 



= 4tt [ 
Jo 



M (r) 
respectively. 



, , 2 dV GM (r)x 

po(r)r dr and - — — = , 

ax r 6 



(D7) 



30 P. O. Vandervoort 



For the radial oscillations of a sphere, it is convenient to construct the Eulerian perturbation of the gravitational accel- 
eration in terms of the Eulerian perturbation of the density pi(x,t) and the hydrodynamical Lagrangian displacement £(x,t) 
(see equations [109]). The perturbation pi is spherically symmetric, so, by Newton's theorem and then Gauss' theorem, the 
Eulerian perturbation of the gravitational acceleration reduces to 



ox r 6 /, , „ 

O \x\ <r 



dx = — 4nr 2 Gpo(r)— ■ £(x, t) = AitOm* f Ax(x,v,t)fo(x,v)dv 



I- 



(D8) 



In the reduction of equation (D8), we have noted that £ lies in the radial direction under conditions of spherical symmetry, 
we have identified i/rasa unit vector in the radial direction, and we have substituted from the second of equations (109) for 

Starting with the definition of the operator Aa{} in equation (3), we insert the Eulerian perturbation of the gravitational 
potential in accordance with equation (29), we make use of equations (D7) and (D8) for — dVo/dx and — dVi/dx, respectively, 
and we find, after additional reductions which are straight-forward, that 

d ( GMo (r) " 



a r t.,i \ 1 GMo(r) . . x 

Aa{z a M{m,n,p)\ — ^-^z^M(m,n,p) — — 



dr y 



/' 



x ■ z a M(m, n,p) + AnGm* / z r7 M(m,n,p)fodv. (D9) 



APPENDIX E: INTEGRALS OVER THE PHASE SPACE 

The applications of the matrix method described in Sections 5 and 6 of this paper are limited to systems in which the 
unperturbed distribution function is a function only of the energy of a star and the velocity distribution is accordingly 
isotropic. Consider in that case the integrals over the velocity space on the right-hand side of equation (D9) and within 
the integrals over the phase space on the right-hand sides of equations (Bl) and (D6). The integrands can be expressed 
as polynomials in the components of v in virtue of equations (7), (33), and (35). Inasmuch as the unperturbed velocity 
distribution is isotropic, integrations over direction in the velocity space are readily performed, and the integrations over the 
velocity space reduce to the evaluation of the even moments 

G(r, l i)=m,J\v\"Mx,v)dv (/i = 0,2,...) (El) 

of the velocity distribution. It is a consequence of the collisionless Boltzmann equation that the moments satisfy the hierarchy 
of ordinary differential equations 

A G(r ^) = _ (/i+ i) G ( r)M _2)^° ( M = 2,4,....). (E2) 

Thus, in the evaluation of the integrals on the right-hand sides of equations (Bl) and (D6), we calculate the moments 
G(r, p,) of the velocity distribution by integrating equations (E2) with the boundary conditions G(r, (j,) — at r — R, where 
R is the radius of the unperturbed system. We note in this connection that the zeroth moment of the velocity distribution is 
G(r,0)=m* / f (x,v)dv = po(r). 

The evaluation of the moments of the velocity distribution reduces the integrands on the right-hand sides of equations 
(Bl) and (D6) to functions of the radial coordinate r = \x\. Accordingly, these integrations over the configuration space reduce 
to quadratures over the radial coordinate r. Those quadratures can be reformulated in terms of the integration of associated 
ordinary differential equations with appropriate boundary conditions. The integration of these equations can be carried out 
simultaneously with the integration of equations (E2) described above. In this work, those integrations have been carried out 
with the aid of a standard Runge-Kutta integrator of the fourth order (Press et al. 1986). 
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